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[i]  We  study  the  tidal  history  of  an  icy  moon,  basing  our  approach  on  a  dissipation  model, 
which  combines  viscoelasticity  with  anelasticity  and  takes  into  account  the  microphysics  of 
attenuation.  We  apply  this  approach  to  Iapetus,  the  most  remote  large  icy  moon  in  the 
Saturnian  system.  Different  authors  provide  very  different  estimates  for  Iapetus’s  despinning 
timescale,  by  several  orders  of  magnitude.  One  reason  for  these  differences  is  the  choice 
of  the  dissipation  model  used  for  computing  the  spin  evolution.  As  laboratory  data  on 
viscoelastic  properties  of  planetary  ices  are  sparse,  many  studies  relied  on  dissipation  models 
that  turned  out  to  be  inconsistent  with  experiment.  A  pure  water  ice  composition,  generally 
assumed  in  the  previous  studies  of  the  kind,  yields  despinning  times  of  the  order  of  3.7  Gyr 
for  most  initial  conditions.  We  demonstrate  that  through  accounting  for  the  complexity 
of  the  material  (like  second-phase  impurities)  one  arrives  at  despinning  times  as  short  as 
0.9  Gyr.  A  more  exact  estimate  will  remain  unavailable  until  we  leam  more  about  the 
influence  of  impurities  on  ice  dissipation.  By  including  the  triaxial-shape-caused  torque, 
we  encounter  a  chaotic  behavior  at  the  final  stage  of  despinning,  with  the  possibility  of 
entrapments  in  the  intermediate  resonances.  The  duration  of  these  entrapments  turns  out  to 
be  sensitive  to  the  dissipation  model.  No  long  entrapments  have  been  found  for  Iapetus 
described  with  our  laboratory-based  dissipation  model. 
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1.  Motivation 

[2]  With  its  albedo  dichotomy,  the  largest  nonhydro  static 
oblateness  observed  in  the  Solar  system,  a  20  km  high  and 
150  km  wide  equatorial  ridge,  an  almost  8°  inclination  to 
Saturn’s  equator  (and,  accordingly,  to  the  orbits  of  other  large 
Saturnian  moons),  Iapetus  is  one  of  the  most  puzzling  objects 
of  the  Solar  System. 

[3]  Although  currently  Iapetus  is  locked  in  a  1 : 1  spin-orbit 
resonance  with  Saturn  and  rotates  at  a  period  of 79.33  days,  its 
original  state  was  most  likely  one  of  more  rapid  rotation,  as 
can  be  deduced  from  Iapetus’  shape.  The  33  km  difference 
between  its  equatorial  and  polar  radii  [ Thomas  et  al. ,  2007] 
has  been  interpreted  as  the  fossil  shape  of  the  satellite  when  its 
spin  period  was  between  15  and  16  h  [ Castillo-Rogez  et  al., 
2007]  Theoretical  models  of  satellite  accretion  (summarized 
by  Castillo  -Rogez  et  al.  [2007])  suggest  that  the  spin  period  at 
the  end  of  accretion  could  be  as  short  as  7  to  10  h. 

[4]  Most  large  moons  of  the  giant  planets  are  believed 
to  have  achieved  the  1:1  spin-orbit  resonance  with  their 
primaries  not  long  after  formation.  The  despinning  time  of  a 
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satellite  scales  as  its  distance  to  the  planet,  to  the  minus  sixth 
power,  and  is  normally  expected  to  lie  within  hundreds  of 
thousands  to  several  millions  of  years. 

[5]  Assuming  a  constant  quality  factor  Q  equal  to  100, 
Peale  [1977]  obtained  for  Iapetus’  a  tidal  despinning  time  of 
order  10  Ga,  much  longer  than  for  other  large  moons.  As 
demonstrated  by  Aleshkina  [2009],  the  use  of  a  dissipation 
model  with  Q  inversely  proportional  to  the  tidal  frequency 
inflates  the  despinning  time  by  two  more  orders  of  magnitude. 
Under  the  initial  conditions  chosen  by  Aleshkina  [2009], 
Iapetus  traverses  the  low-order  resonances  nonstop.  At  the 
same  time,  it  was  acknowledged  by  Aleshkina  [2009]  that 
different  initial  conditions  might  cause  a  temporary  capture  in 
some  of  the  resonances,  in  which  case  the  despinning  time- 
scale  would  be  even  longer. 

[6]  The  construction  of  more  adequate  models  of  tidal 
dissipation  in  Iapetus  and  other  satellites  is  thus  needed. 
First  and  foremost,  these  models  should  include  the  actual 
dependence  of  the  damping  rate  upon  the  tidal  frequency.  As 
demonstrated  by  Efroimsky  and  Lainey  [2007],  employment 
of  a  realistic  dissipation  law  based  on  experimental  mea¬ 
surements  can  considerably  influence  the  timescale  of 
dynamical  evolution  of  the  Solar  system  objects  and  is  also 
relevant  to  exoplanetary  systems.  At  a  more  comprehensive 
level,  a  model  should  also  account  for  the  feedback  between 
the  thermodynamical  state  and  the  forcing  stress  and  fre¬ 
quency.  For  example,  many  minerals  (including  ices)  are 
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known  to  become  more  dissipative  with  increasing  temper¬ 
ature  (i.e.,  with  decreasing  viscosity). 

2.  Statement  of  Purpose 


each  tidal  change  of  the  potential,  Wj,  entails  a  linearly  pro¬ 
portional  deformation  of  the  primary’s  shape.  The  deforma¬ 
tion  will,  in  its  turn,  amend  the  potential  of  the  primary  with 
a  linear  adjustment  Up. 


[7]  Behavior  of  icy  moons  under  high  tidal  stress,  such  as 
Europa  or  Enceladus  (where  the  tidal  stress  is  of  the  order  of 
105  Pa),  has  received  much  attention,  since  for  these  bodies 
tides  are  easily  identified  as  a  major  driver  of  endogenic 
activity.  At  the  same  time,  tidal  dissipation  and  its  impact  on 
the  internal  evolution  of  objects  subject  to  very  low  stressing 
still  await  exploration.  This  is  why  we  address  Iapetus,  an  icy 
moon  experiencing  tidal  stresses  of  the  order  of  ICE  Pa. 

[8]  We  construct  a  new  dissipation  model  which  rests 
on  our  laboratory-based  understanding  of  the  viscoelastic 
properties  of  ice.  While  no  attenuation  measurements  on  icy 
materials  have  yet  been  obtained  under  conditions  exactly 
identical  to  those  on  Iapetus,  our  approach  is  based  on  a  large 
bulk  of  relevant  experimental  observations  and  theoretical 
considerations  available  by  now.  The  approach  includes 
quantifying  the  friction  rate  as  a  function  of  the  deformation 
mechanisms  expected  to  act  in  Iapetus  under  the  assumption 
of  low  porosity.  In  our  computations,  we  trace  the  time 
evolution  of  the  temperature  and  the  excitation  frequency, 
and  thereby  of  the  viscoelastic  properties  of  the  material.  Our 
simulations  are  based  on  the  despinning  theory  developed  by 
Efroimsky  and  Williams  [2009],  combined  with  the  geo¬ 
physical  model  by  Castillo -Rogez  et  al.  [2007],  Specifically, 
we  amend  that  model  with  a  realistic  dissipation  law  which 
relies  on  experimental  measurements.  We  also  consider  the 
effect  of  Iapetus’  triaxial  shape  (triaxiality)  on  despinning. 
Finally,  we  explore  the  important  question  of  the  evolution  of 
the  tidal  torque  on  approach  to  a  spin-orbit  resonance. 

[9]  In  section  3,  we  recall  the  meaning  of  linearity  and 
discuss  the  frequency  dependencies  of  the  quality  factor,  tidal 
Love  numbers,  and  phase  lag.  In  section  4,  we  present  a 
fonnula  for  the  tidal-despinning  rate.  Section  5  summarizes 
the  previously  performed  research  on  Iapetus’  tidal  despin¬ 
ning.  Section  6  discusses  the  rationale  for  introducing  a  new 
dissipation  model  consistent  with  the  temperature  and  stress 
conditions  expected  in  Iapetus,  supported  by  geophysical 
considerations.  In  section  7,  we  assemble  our  model  in  its 
entirety  and  present  the  results  of  numerical  calculations.  We 
also  probe  the  possible  role  of  the  triaxiality.  Conclusions  are 
drawn  in  section  8.  In  Appendix  A,  we  explain  why  the 
realistic  dissipation  models  (with  Q  scaling  as  a  positive 
power  of  frequency)  do  not  entail  infinities  in  the  expressions 
for  the  tidal  torque.  We  thus  refute  a  popular  fallacy  that  such 
rheologies  yield  diverging  torques  in  the  zero-frequency 
limit. 

3.  The  Quality  Factor,  Phase  Lag,  and  Love 
Numbers 

3.1.  Tidal  Deformation  of  a  Homogeneous  Spherical 
Primary 

[10]  A  stationary  potential  W(R,  R*),  generated  by  a  sec¬ 
ondary  residing  at  R *  outside  the  primary,  can  be  expanded, 
at  each  point  R  on  the  primary’s  surface,  over  the  Legendre 
polynomials  P/(cos7).  Here  7  denotes  the  angular  separation 
between  R*  and  R,  both  vectors  coming  out  of  the  primary’s 
center.  Stationary  deformation  is  linear  if  in  the  expansion 


Ui(R)  ~  Wi(R,  R*).  (1) 


[11]  For  a  spherical  primary  with  mean  equatorial  radius  R, 
the  potential  at  degree  /  decreases  outside  the  surface  as 
1  lr>+  '.  Hence,  the  linearity  assertion  results  in 

/  d\  J+l 

t/;(R)  =£;(-)  R)(R.R*),  (2) 


kj  being  the  Love  number,  R*  =  (/■*,  (j>*,  A*)  being  the  radius, 
latitude,  and  longitude  of  the  tide-raising  secondary,  r  = 
(. R ,  </>,  A)  being  the  coordinates  of  a  surface  point,  and  r  = 
(r,  <f>,  A)  being  those  of  an  exterior  point  located  above  it  at 
a  radius  r>R. 

[12]  Leaving  consideration  of  radially  stratified  objects  for 
section  3.5,  we  begin  with  the  simpler  case  of  homogeneous 
bodies.  For  a  homogeneous  incompressible  spherical  primary 
of  density  p,  surface  gravity  g,  radius  R,  and  static  rigidity  p 
(or  the  static  compliance  J=  Up),  the  static  quadrupole  Love 
number  is  given  by  [ MacDonald ,  1964] 


1 

1+^2 


where 


51  J-' 

8  7T  7  p2  R2  ' 


k 


_  19 

Al  ~  2  pgR 


57  p 

8  7T  7  p1  R2 


(3) 


7  =  6.7  x  10  1 1  m3  kg  1  s  being  Newton’s  gravity  constant. 
For  an  arbitrary  /,  the  general  formula  is: 


k,  = 


2(/  —  1)  1  +A,' 


{2k  -j-  4/  3 )/i 

where  At  = - — 

IgpR 


3(2/2  +  4/  +  3)/i 
4  /  n  G  p2  R2 


(4) 


[13]  All  the  above  pertains  to  stationary  loads  and, 
accordingly,  static  Love  numbers.  A  realistic  tide  varies  in 
time  and  can  be  expanded  into  modes  given  by  expression 
(37)  below.  Historically,  the  modes  Lohnpq  are  numbered  with 
four  integers  [Kaida,  1964]:  Impq,  where  /  >  2.  The  modes’ 
absolute  values,  Ximpq  =  \^impq\,  are  the  actual  physical 
frequencies  of  stresses  and,  consequently,  of  strains.  Due 
to  internal  friction,  a  spectral  component  of  the  strain  lags 
behind  the  appropriate  component  of  the  stress  by  a 
frequency-dependent  phase  shift  €impq(x),  the  functional 
form  of  the  frequency  dependence  being  different  for  differ¬ 
ent  values  of  /  but  independent  of  m,  p,  q.  This  is  why  the 
customary  notation  €lmpq  may  be  substituted  with  a  nota¬ 
tion  depicting  the  situation  more  accurately:  efycimpq)-  The 
Love  numbers  kj(xi„lpq),  too,  depend  upon  the  frequency,  the 
dependencies  looking  different  for  different  Is  and  bearing  no 
dependence  upon  m,p,  q.  (The  situation  changes  in  bodies  of 
a  triaxial  shape.  There,  coupling  between  spherical  hannonics 
renders  the  Love  numbers  and  lags  whose  expressions  via  the 
frequency  depend  on  m,p,  q  [Dehant,  1987a,  1987b;  Smith, 
1974].) 
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[14]  To  keep  the  model  linear  under  varying  load,  we 
assume  that  for  deformation  at  frequency  Ximpq,  the  appro¬ 
priate  Love  numbers  kfximpq)  and  lags  tfximpq )  depend  on 
Ximpq ,  and  not  on  the  other  frequencies  in  the  spectrum,  nor 
on  the  deformation  magnitude  at  this  or  other  frequencies. 

[15]  When  lagging  is  linear,  the  average  dissipation  rate 
(E(x))  and  a  one-cycle  energy  loss  A Ecvcle(x)  at  a  frequency 
X  can  be  expressed  with  the  empirical  relations 

(£(X))  =  -X^g^  and  A£0,cte(X)  =  -2.^M. 

(5) 

If  Epeak{x)  is  the  peak  energy  stored  at  the  frequency  x,  the  Q 
factor  is  related  to  the  lag  via 


the  complex  amplitudes  being: 

0&{x)  =  ^o(x)  d'p"{x\  Mx)  =  Mx)  e^u(x).  (10) 

[22]  At  each  frequency  x,  the  initial  phases  <fi„(x)  and 
tpu(x)  can  be  picked  so  that  the  amplitudes  ay(x)  and  Uy(x) 
are  nonnegative. 

[23]  In  the  case  of  a  continuous  spectrum,  the  sums  will 
become  integrals: 

PCO  /*  GO 

f 7<p{t )  =  /  ay(x)  e,x‘  d\  and  uy(t)  =  /  iiy(x)  e‘x‘  d\, 

Jo  Jo 

(11) 


Q  1  =  sin|e|  (6) 

and  not  Q  1  =  tan  |  e  as  often  presumed.  Indeed,  the  latter 
definition  is  widely  used  in  the  experimental  literature  but  is 
not  directly  applicable  in  the  present  case.  If  Epeak(x)  is 
defined  as  the  peak  work,  the  corresponding  Q  factor  will  be 
related  to  the  lag  in  a  more  complicated  manner,  as  demon¬ 
strated  by  Efroimsky  and  Williams  [2009]: 


tan|e| 

1  _  N)tanH 


(7) 


[16]  Efroimsky  and  Williams  [2009]  were  inaccurate  when 
they  called  Epeak(;x)  the  peak  energy.  However,  their  calcu¬ 
lation  of  Q  was  carried  out  with  the  understanding  that 
Epealfx )  is  the  peak  work. 

[17]  In  the  limit  of  small  e,  expression  (7)  becomes 

Q1  =  sin|e|  +  0(e2),  (8) 

so  definition  (7)  makes  1  IQ  a  good  approximation  to  sin  e  for 
small  lags  only. 

[is]  Both  definitions  also  render  Q  — *  0  for  e  — *  7r/2. 

[19]  The  Darwin-Kaula  expansion  of  tides  contains  not 
the  inverse  quality  factors,  but  sines  of  the  phase  lags  (see,  for 
example,  fonnula  (102)  of  Efroimsky’  and  Williams  [2009]  or 
our  formula  (35)).  Thus,  the  despinning  theory  used  in  this 
study  will  operate  with  sin  e  and  not  with  1  IQ,  as  the  defi¬ 
nition  of  Q  is  genuinely  ambiguous. 

[20]  Getting  back  to  the  Love  numbers,  we  would  recall 
that  a  potential  proportional  to  Pj( cosy)  must  be  decreasing 
outside  the  primary  as  r  (l  +  1  \  Then  we  see  that  the  tidal 
dynamics  is  determined  almost  exclusively  by  the  principal 
Love  number  k2. 


where  we  omit  the  symbol  7 Ze  and,  for  evident  physical 
reasons,  integrate  over  positive  x  only.  Whenever  necessary, 
the  frequency  will  be  assumed  to  approach  the  real  axis  from 
below:  lm(x)  <  0,  so  (11)  will  be  a  Laplace  rather  than  a 
Fourier  transform. 

[24]  Under  the  assumption  of  incompressibility,  the  volu¬ 
metric  part  of  the  strain  is  nil.  Hence  in  what  follows  we  shall 
deal  only  with  the  deviatoric  components.  In  a  linear  material 
with  memory,  the  strain  is  rendered  by  the  compliance 
operator  J(t): 


2  uy(f)  =  j(t)  ay  =  /  J(t)  & y(t  —  t)  dr 
Jo 

=  f  J(t  —  t')  dt\ 

J  —OO 


(12a) 


the  overdot  denoting  time  derivative.  The  kernel  J(t  -  t'), 
termed  compliance  function  or  creep  response  function 
[Karato,  2008],  describes  the  mechanical  behavior  of  a 
medium  (i.e.,  its  capacity  to  comply  )  in  response  to  an  applied 
stress,  and  thus  characterizes  its  memory.  A  stress  increment 
da^it')  =  a^ft')dt'  at  the  time  t'  contributes  to  the  strain  at  a 
later  time  as  follows:  2  du^ft)  =  J(t  -  t')da^v(t'). 

[25]  After  integration  by  parts,  the  operator  acquires  the 
form  of 


2  Uy(t)  =  J(t)  ay  =  J( 0)  a<y{t)  -  j(oo)  ay(-co) 


+ 


j(t  —  t')  ay(t')  dt' 


(12b) 


2  Uy(t)  =  -J(a o)  ay(-co) 

[‘  d 

+  /  —[J{t-t')-J(0)+J(0)Q(t  —  t’)\ay(t')dt', 

(12c) 


3.2.  Complex  Notations  for  the  Stress,  Strain,  Rigidity, 
and  Compliance 

[21]  Insofar  as  the  linear  approximation  remains  valid,  the 
stress  tensor  ay  and  the  strain  tensor  Uy  can  be  expanded 
into  Fourier  series 

/(0  =  TZe[ay(x)e,x‘] ,  uy(t)  =  ^  TZe[uy(x)e'x'} , 

X  X 

(9) 


0(7  -  t')  denoting  the  Heaviside  step  function,  which  van¬ 
ishes  for  t  ~  t'  <  0  and  assumes  the  value  of  unity  for  t  —  t'  >  0. 

[26]  Assuming  that  the  stress  ay{~ 00)  infinitely  long  ago 
was  zero,  we  drop  the  term  containing  the  relaxed  compliance 
J( 00).  However,  the  unrelaxed  compliance,  J( 0),  is  always 
important,  as  it  describes  an  immediate  response  to  forcing.  It 
can,  though,  be  incorporated  into  the  kernel.  As  we  just  saw 
in  (12c),  this  can  be  done  if  we  accept  that  the  unrelaxed 
compliance  enters  the  compliance  function  not  as  J(t  -  t')  = 
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J(0)  +  . . .  but  as  J{t  ~  t')  =  ~t')+  ....  Hence,  in  what 

follows  we  shall  always  write 


2  =  /  j(t  —  t ')  dt'. 


(13) 


implying  that  J(t  -  t')  includes  J( 0)  0(7  -  t ')  instead  of  J(0). 

[27]  Complementary  to  the  shear  compliance  operator  J(t) 
is  the  operator  of  shear  stress  relaxation  which  may  as 
well  be  called  the  operator  of  rigidity: 


pco 

cr$ „(f)  =  2  ft(t)  uc/j  =  2  h(t)  iifyit  —  r)  dr 

Jo 

=  f  2  n(t  —  t ')  iicyi?)  dt'.  (14) 

J  —00 


Its  kernel,  n(t  -  t'),  named  the  stress  relaxation  function,  may 
include  the  term  /i(0)O(t  - 1')  responsible  for  the  elastic  stress. 
Here  /i(0)  is  the  unrelaxed  rigidity,  an  inverse  to  the  unrelaxed 
compliance:  /i(0)  =  1/J(0).  The  kernel  may  also  include  the 
viscous  term  2 r]6(t  -  t'),  with  8(t  ~  t')  denoting  the  delta 
function,  and  >7  being  the  viscosity. 

[28]  In  the  presence  of  viscosity,  ft  thus  becomes  an  inte- 
gro-differential  operator,  and  not  an  integral  operator  like  J . 
This  makes  the  integration  by  parts  like  (12b)  and  (12c) 
impossible.  Accordingly,  we  cannot  write  the  stress  as 
aCv{t)  =  2/(_00  (i(t  ~  t')u0it')dt’. 

[29]  Introducing  the  complex  compliance  J(x)  through 


j: 


J(x)e’XTdx  =  J{t)  where  J(x)  =  /  j(r)e  ‘XTdr, 

Jo 

(15) 


3.3.  The  Love  Operators  and  the  Complex  Love 
Numbers 

[30]  Consider  a  homogeneous  incompressible  primary 
subject  to  a  stationary  disturbance  by  a  secondary.  The  /th 
spherical  harmonic  of  the  disturbance-caused  increment  of 
the  primary’s  exterior  potential  is  traditionally  denoted  as 
U/(R).  This  increment  is  related  to  the  /th  spherical  harmonic 
W](R,  R)  of  the  perturbing  exterior  potential  via  relation  (2). 

[31]  For  time-dependent  deformations,  the  Love  numbers 
become  integral  operators,  so  relation  (2)  must  be  generalized 
to 


i+\ 


I/,(R,0=(-1  ki(t)  Wi(R.  R*.  t'). 


(21) 


Under  weak  deformations,  it  is  natural  to  assume  this  operator 
to  be  linear,  so 


O/(R.0=  (?) 


'R\ 


/+1  /»a 

Jo 

/+!  pt 


k/(r)JV/(R,  R*,  t  —  r)dr 
ki(t-t')Wi(R.R*,t')dt'. 


(22a) 


Integration  by  parts  provides  an  equivalent  expression: 

/  d\  t+1  rco 

C/(R,t)=f-J  J  ki{r)Wi{R,R*,t-T)dT 


R\ 

r 


^+1  pt 


> 


ki(t-t')Wi(R,R*,t')dt'. 


(22b) 


and  defining  the  complex  rigidity  /z(x)  as 

P(x)  =  1  (16) 

one  can  interconnect  the  complex  amplitudes  (11)  through 

0fr(x)  =  2  Hx)  «^(x),  2  uc„(x)  =  J(x)  ff^(x)-  (17) 

(The  presence  of  the  viscous  term  2 r/8(t  -  t')  in  the  kernel 
fi(t  ~  /'),  prevents  us  from  defining  j2(x)  via 

00  poo 

Ji(x)e'XTdx  =  where  Ji(x)  =  /  /i(r)e‘‘XT6?r, 

-00  Jo 


Just  as  in  expressions  (12b)-(12c)  for  the  compliance  oper¬ 
ator,  in  (22b)  we  obtain  the  boundary  terms  ki(0)W(t)  and 
— k/(co)W(—  00).  While  setting  W(—  co)  =  0  justifies  our  neglect 
of  the  latter  term,  the  former  term  may  be  absorbed  into  the 
kernel  by  incorporating  k/(0)O(t  -  t')  into  kt(t  -  t'). 

[32]  Expressions  (22)  coincide  with  (2)  in  the  limit  of  a 
perfectly  elastic  body.  Indeed,  in  this  case  kt(t  -  t')  contains 
nothing  but  the  immediate  reaction  term  A'/(O)0(t  -  t'),  while 
the  derivative  of  k 7  becomes:  ki(t  -  t ')  =  k/(0)8(t  -  t').  Some¬ 
times  the  time  derivatives  At/(t)  are  called  Love  functions,  a 
term  coined  by  Churkin  [1998]. 

[33]  In  terms  of  the  complex  spectral  components  of  U i('X) 
and  fV/(x),  one  can  rewrite  (22)  as 


lest  we  have  to  deal  with  a  derivative  of  the  delta  function.) 
Writing  the  complex  rigidity  and  compliance  as 


U/(x)  =  k/(x)  W/(x),  (23) 


A»(x)  =  \d(x)\  exp[«5(x)]  and  J(x)  =  |/(x)|  exp[-5(x)], 

(18) 


one  obtains 


k/(x)  being  related  to  Ay(r)  =  dki(r)/dT  via 


A/(r)  =  /  ki(x)  e‘XTdx 


(24) 


tan  S(x) 


lm[J(x)\  _  Tm\p{xj] 
ne[j{x)}  Mrtx)] ' 


(19) 


Evidently,  the  phase  angle  8(x)  is  a  measure  of  lagging  of  the 
strain  at  the  frequency  x  behind  the  stress  at  this  frequency: 


VN(x)  =  w(x)  -  S{x)-  (20) 


(the  idea  of  making  the  Love  number  complex  was  pioneered, 
likely,  by  Munk  and  MacDonald  [1960]  and  Zschau  [1978]). 
Insofar  as  alternating  deformation  remains  linear,  it  obeys  the 
elastic -viscoelastic  analogy,  a  correspondence  principle  that 
relates  a  solution  of  a  linear  viscoelastic  boundary  value 
problem  to  an  analogous  problem  of  elastic  body  mechanics, 
with  the  same  initial  and  boundary  conditions  [Haddad, 
1995].  As  a  result  of  this  correspondence,  the  algebraic 
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equations  for  the  Fourier  (or  Laplace)  components  of  the 
strain  and  stress  in  the  viscoelastic  problem  mimic  the 
equations  interconnecting  the  strain  and  stress  in  the  appro¬ 
priate  elastic  problem.  On  these  grounds,  the  operational 
moduli  can  be  manipulated  algebraically  in  the  same  way  as 
their  elastic  counterparts.  For  this  reason,  the  complex  Love 
numbers  are  linked  to  the  complex  rigidity  in  the  same 
manner  as  the  static  Love  numbers  depend  upon  the  static 
rigidity.  This  correspondence,  however,  remains  in  force 
insofar  as  the  linearity  condition  is  met.  As  explained  in 
section  6,  lapetus  meets  this  requirement:  the  tidal  stress  in  it 
is  low  enough  to  keep  deformation  linear  (the  dissipation  thus 
being  dominated  by  diffusion  of  defects  through  the  lattice,  a 
linear  mechanism).  At  the  same  time,  the  developed  theory 
may  be  too  crude  for  the  icy  satellites  subject  to  tidal  stress  of 
magnitudes  greater  than  ~0.1  MPa  (e.g.,  Europa),  inside 
which  the  nonlinear  mechanisms  of  grain  boundary  sliding  or 
of  dislocation  creep  are  likely  to  be  dominant.  The  applica¬ 
bility  ofthe  correspondence  principle  to  those  objects  remains 
questionable. 

[34]  According  to  the  correspondence  principle,  the  alge¬ 
braic  dependence  of  a  complex  Love  number  on  the  complex 
rigidity  mimics  the  relation  between  their  static  counterparts. 
The  quadrupole  Love  number  of  a  homogeneous  near- 
spherical  body  will  look: 


where  e/(x)  is  the  phase  lag  at  the  tidal  frequency  It  is 
defined  through 


tan  e;(x) 


Zm[k,{x)\ 

Ke[h{x)\ 


(29) 


or,  equivalently,  through 

\h(x) I  sin ei(x)  = -lm[k,(x)\ 

3 _ -A,Jlmp(X)] 

2(1  -  1)  (7 Ze[J(X)\  +  A,  jf+(Xm[J(X)]f  ' 

(30) 

The  importance  of  this  formula  stems  from  the  fact  that  the 
product  kt  sin  e7  =  |Ay(x)|sin  e/(x)  emerges  in  the  Darwin- 
Kaula  expansion  for  the  tidal  potential.  Hence  it  is  this 
product  (and  not  kj/Q)  that  enters  the  tidal  force,  torque,  and 
dissipation  rate.  This  circumstance  paves  the  way  from  a 
dissipation  model  to  the  despinning  history.  Indeed,  the  dis¬ 
sipation  model  implies  a  certain  dependency  of  /  on  the  fre¬ 
quency  X-  The  form  of  this  dependency  defines  the  frequency 
dependence  of  the  Love  numbers,  kj (x).  The  latter  define,  for 
each  1,  the  functional  form  of  the  dependencies  (30)  which, 
in  their  turn,  are  directly  employed  in  the  tidal  theory. 


ki{x) 


3  1  _  3  1  _  3  i 

2  j  19/j(x')  2  1  19  p  Ji(x)  2  1  +  A2  p(x)// 

2pgR  2pgR  p 

(25a) 


or,  in  terms  of  the  complex  compliance: 


ki(x) 


3  1 

2  1  +A2  J fJ (x) 


3  _  /(x) 

2  J(x)  +  A2J 


(25b) 


where  the  dimensionless  factor  is  given  by 

19/u  57/t  _  57  J~l 

2  pgR  8  7r  7  p2  R2  8  7r  7  p2  R2’ 


(26) 


7  being  Newton’s  gravity  constant.  Formally,  this  expression 
for  A2  coincides  with  its  static  counterpart  (3).  A  difference, 
however,  exists.  While  in  (3)  letters  p  and  /  signified  the  static 
rigidity  and  compliance,  in  (26)  they  may  stand  for  any 
benchmark  values  obeying  p  =  1  //.  The  reason  for  this  is  that 
the  product  A2J  entering  (25)  does  not  in  fact  depend  upon  / 
or  p.  It  may  be  reasonable  to  use  the  unrelaxed  values  p  =  p( 0) 
and  J  =  J( 0).  However,  this  choice  will  not  be  obligatory, 
because  in  some  rheological  models  the  unrelaxed  moduli 
may  be  zero  or  infinite. 

[35]  For  an  arbitrary  /,  (25)  will  become: 

3  1  3  J(x) 

’  2(/-l)  l+Ai  J/J(X)  2(I-l)J(x)+A,J’  {  ’ 


with  the  above  caveat  concerning  the  expression  of  A/  via 
p  or  J. 

[36]  The  complex  Love  number  can  be  naturally  written  as 


ki(x)  =  ne\k,{x)\  +  i1m[ki(x)\  =  \ki(x)\  e  K,(x)  (28) 


3.4.  The  Case  of  an  Incompressible  Homogeneous  Body 

[37]  lapetus’  rigidity,  being  temperature-  and  porosity- 
dependent,  is  likely  to  have  varied,  over  the  satellite’s  history, 
within  the  range  of  p  ~  4  x  108  to  4  x  109  Pa.  If  we  model 
lapetus  with  an  incompressible  and  homogeneous  sphere  of 
density  p  ~  103  kg  nfr3  and  radius  R  ~  7.5  x  105  m,  these 
numbers  put  the  value  of  A /  within  the  interval  of  25-250. 
This  tells  us  that 


An 


J 


\J(X)\ 

so  we  can  approximate  (27)  with 

3  J(X)  „  3  J(X) 


>  L 


Hx)  = 


2  J(X)  +  A,J  2  A,J 


o(|7/(/,/)|2), 


(31) 

(32) 


except  in  the  closest  vicinity  of  the  resonance,  where  the  tidal 
frequency  \  approaches  zero,  and  /  diverges  for  some 
mechanical  models,  like,  for  example,  for  those  of  Maxwell 
or  Andrade.  (Recall  that  according  to  (26),  the  dimensionless 
factor  A 1  is  inversely  proportional  to/,  so  the  product ///does 
not  in  fact  depend  upon  /.) 

[38]  Whenever  the  approximate  formula  (32)  is  applicable, 
we  can  rewrite  (29)  as 


tan  e/(x) 


lm\ki(x)\ 

Ke[ki(x)\ 


lm[J(x)\ 

Ke(J(x)\ 


tan£(x),  (33) 


whence  the  phase  lag  of  a  tidal  mode  at  frequency  x  coincides 
with  the  phase  lag  in  a  sample: 

e/(x)  ~  <5(x)  for  ah  l  (34) 

provided  x  is  not  too  close  to  zero  (i.e.,  we  are  not  too  close  to 
a  resonance).  Stated  alternatively,  for  \  not  approaching  zero 
too  closely,  the  component  Ui(x)  of  the  primary’s  potential 
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variation  lags  behind  the  component  W /(x)  of  the  perturbed 
potential  by  the  same  phase  angle  as  the  strain  lags  behind  the 
stress  at  frequency  x  hi  a  sample  of  the  material. 

[39]  Just  as  (30),  so  formulae  (32)  and  (33)  should  be 
employed  with  an  important  caveat  in  mind.  In  reality,  the 
potential  U  and  therefore  also  k  are  functions  not  of  the 
positively  defined  physical  frequencies  x  but  of  the  tidal 
modes  u>,  which  can  be  positive  or  negative.  Our  illegitimate 
use  of  x  instead  of  u>  should  be  compensated  through 
multiplying  the  lag  f-i('Ximpq\  “by  hand,”  with  sgn wlmpq.  So, 
for  example,  (34)  should  actually  look 

e/(x)  ~  5(x)sgnw. 

Leaving  the  (somewhat  tedious)  mathematical  elucidation  of 
this  fact  to  another  paper,  we  would  mention  that  from  the 
viewpoint  of  physics,  the  need  for  the  sign  factor  is  clear. 
Indeed,  both  the  lag  elmpq  and  the  appropriate  term  of  the 
torque  should  change  their  sign  on  crossing  the  lmpq  com- 
mensurability  (while  the  lag  6  in  the  material  should  always 
stay  positive,  as  the  strain  always  falls  behind  the  stress 
causing  it). 

[40]  Condition  (31)  is  obeyed  insofar  as  changes  of  shape 
are  determined  solely  by  the  material  properties,  and  not 
by  self-gravitation  of  the  object.  This  condition  therefore 
restricts  the  applicability  realm  of  (34).  An  even  tougher 
limitation  stems  from  the  fact  that  the  above  treatment  relied 
on  the  homogeneity  assumption.  Realistic  models  of  larger 
bodies  must  account  for  the  radial  variations  in  the  material 
properties  due  to  temperature  gradients  and  possible  com¬ 
positional  stratification.  (Such  stratification  may  emerge  if, 
for  example,  an  icy  moon  has  undergone  internal  melting, 
with  the  ensuing  chemical  differentiation  of  the  interior.  )  For 
lapetus,  thermal  models  suggest  that  compositional  differ¬ 
entiation  was  limited  [Castillo -Rogez  et  al.,  2007].  Still,  the 
impact  of  thermal  variations  on  the  mechanical  properties  of 
the  material  was  substantial  enough.  So  we  cannot  directly 
employ  (34)  but  instead  have  to  derive  the  overall  lag  e 
through  integration  over  layers. 

3.5.  The  Case  of  a  Radially  Heterogeneous  Body 

[41]  The  dissipative  properties  of  layered  objects  are 
explored  by  numerical  integration  of  the  deformation  prop¬ 
agation  from  the  interior  to  the  surface,  over  the  layers  of  the 
body.  In  our  computations,  we  relied  on  the  model  developed 
by  Takeuchi  and  Saito  [1972]  for  an  incompressible  near- 
spherical  object.  The  incompressibility  assumption  is  gener¬ 
ally  adopted  in  the  literature  [e.g.,  Tobie  et  al. ,  2005]  because 
compression  plays  only  a  minor  role  in  the  overall  tidal 
response  of  a  body. 

[42]  The  overall  phase  lag  e  is  computed  via  integration  of  6 
over  the  layers.  Hence,  e  becomes  a  function  of  the  compo¬ 
sition  stratification,  the  physical  structure  (e.g.,  the  porosity 
and  the  presence  of  a  liquid  layer),  and  the  temperature  dis¬ 
tribution.  Computation  of  e  includes  computation  of  the 
profile  of  stresses  and  strains  over  the  layers,  the  compliance 
properties  of  the  medium  in  each  layer  being  different  (though 
always  remaining  linear  (see  section  6)). 

[43]  The  global  deformation  is  inferred  from  the  equation 
of  motion,  as  explained  by  Takeuchi  and  Saito  [1972],  His¬ 
torically,  their  method  was  developed  to  model  a  static  elastic 


deformation  of  a  multilayered  body,  so  the  radial  functions 
and  equation  of  motion  of  Takeuchi  and  Saito  [1972] 
contained  the  static  rigidity.  Thanks  to  the  correspondence 
principle,  the  method  can  be  applied  to  varying  loads  too.  The 
radial  functions  and  equations  of  motion  (as  well  as  the 
constitutive  equations)  are  then  written  in  the  frequency 
domain  and  thus  contain  either  the  complex  rigidity  /7(x)  or 
its  inverse,  the  complex  compliance  J(x).  The  method 
remains  agnostic  of  the  viscoelastic  model  employed,  insofar 
as  the  model  remains  linear  so  the  correspondence  principle 
work. 

[44]  The  complex  tidal  Love  numbers  ki(\)  are  determined 
from  the  radial  functions  integrated  to  the  surface.  Assuming 
that  the  body  is  spherically  symmetric  and  isotropic  (and  that 
neither  of  these  two  assumptions  is  violated  too  badly  by 
convection),  the  deformation  vector  u  generated  by  the  tidal 
potential  can  be  expanded  over  the  spherical  functions  Y\n 
multiplied  by  some  functions  of  radius  (the  latter  functions 
being  different  for  different  layers).  Details  of  numerical 
solution  of  the  resulting  system  of  differential  equations  are 
given  by  Tobie  et  al.  [2005]. 

[45]  The  difference  between  the  thus  obtained  Love  num¬ 
ber  for  an  idealized  spherical  object  and  the  appropriate  Love 
number  for  an  actual,  slightly  triaxial  object  is  of  the  order  of 
the  flattening  [Dehant,  1987a,  1987b;  Smith ,  1974].  Taken 
the  uncertainty  in  our  knowledge  of  the  physical  parameters,  a 
small  nonsphericity  can  then  be  neglected. 

4.  Parameterization  of  the  Despinning  Rate 

[46]  While  the  MacDonald  [1964]  theory  tacitly  sets  the 
quality  factor  inversely  proportional  to  the  frequency,  the 
theory  of  Darwin  [1879]  and  Kaula  [1964]  is  general  enough 
to  accommodate  an  arbitrary  dissipation  model.  In  the  fol¬ 
lowing,  we  shall  employ  their  theory  and  shall  treat  lapetus  as 
the  primary  and  Saturn  as  the  secondary.  Our  first  step  will  be 
to  bring  in  the  conventional  expression  for  despinning  rate. 
As  it  is  common  to  write  this  expression  with  the  Q  factor  in 
the  denominator,  our  second  step  will  be  to  recall  that  in 
reality,  1/(7  is  but  an  approximation  to  the  correct  factor  sin  e, 
as  can  be  easily  seen  from  equation  (35).  For  sufficiently  large 
Q,  the  approximation  works.  However,  experimental  data 
show  that  under  the  conditions  expected  in  icy  moons,  the  Q 
factor  of  ice  can  become  close  to  and  even  lower  than  unity 
[e.g.,  McCarthy  et  al.,  2008].  Thus,  while  the  approximation 
is  justified  for  telluric  bodies,  it  is  not,  generally,  applicable  to 
icy  moons.  So  we  shall  have  to  stick  to  the  phase  lag  e. 

[47]  As  the  inverse  quality  factor  is  (up  to  27r)  the  relative 
energy  loss,  one  may  enquire  if  the  values  of  Q  can  at  all  be 
close  to  or  below  unity.  For  seismic  waves  (which  are  similar 
to  an  underdamped  undriven  oscillator),  the  answer  is  nega¬ 
tive.  It  is,  however,  affirmative  for  tides,  for  these  are  similar 
to  a  driven  (and,  in  some  cases,  overdamped)  oscillator.  The 
quality  factor  falling  short  of  unity  indicates  that  the  eigen- 
frequencies  get  damped  away  during  less  than  one  vibration 
and  that  the  motion  continues  only  due  to  the  driving  force. 

[48]  Our  third  step  will  be  to  recall  that  e  is  the  lag  between 
the  reaction  of  the  primary’s  shape  and  the  external  load.  It  is, 
however,  the  phase  lag  8,  whose  frequency  dependence  is 
measured  in  the  lab.  As  the  body  is  multilayer,  we  cannot 
approximate  e  with  5  but  will  have  to  average  over  layers  to 
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obtain  the  overall  -Tm[k/(x)\  =  £/(x)sin  £/(x)  for  the  body  as  a 
whole. 

4.1.  Despinning  Torques  in  the  Darwin-Kaula  Theory 
of  Bodily  Tides 

[49]  In  the  presence  of  only  one  secondary  (here, 
Saturn)  located  at  the  position  r  =  (r,  <fi ,  A)  relative  to 
the  primary  (here,  Iapetus),  the  tidal  torque  acting  on  the 
primary  reads 


^  =  E2  yM^'+'a-2'-2 
1=2 


(/  —  m)\ 


m 


E  Flp(i) 

P= 0 


00 

■  E  Glq(e)klsinelmpq  +  T 

q=— 00 


00  / 

=  -E2^2,+1«-2,-2E 

1—2  m= 0 


(/  —  m)\ 

-  YYl 

(l  +  m)\ 


p=0 


00 

■  E  G^(e)T^'(e*«)]  +  ?,  (35) 

q=—co 


must  be  derived  from  the  expression  for  the  complex  Love 
number  J  fximpq)-,  which  in  its  turn  is  deduced  from  the  fre¬ 
quency  dependence  J(ximPq)-  In  the  case  of  a  homogeneous 
body,  the  latter  is  rendered,  for  /  =  2,  by  (25)  and  is  often  well 
approximated  by  (32). 

4.2.  The  Formula  for  Despinning  Rate 
and  Its  Limitations 

[51]  As  evident  from  expansion  (35),  the  contribution 
of  tidal  mode  uJi,npq  to  despinning  is  proportional  to  kjfx) 
sin  eimpq(x),  where  x  =  Ximpq  =  Wi,„Pq\-  The  principal  tidal 
mode 


W2200  =  2  («  —  9)  (39) 

generates  the  physical  flexure  frequency 

X2200  =  2|ft  —  d\  =  2(8  —  ft)  (40) 

and  the  phase  lag 


M2ec  being  the  mass  of  the  secondary  (Saturn),  the  sum 
standing  for  the  constant  (independent  from  the  mean 
anomaly  A4)  part  of  the  torque,  M^ec  being  the  mass  of  the 
tide-raising  secondary  (whose  role  in  this  setting  is  played 
by  Saturn),  and  T  denoting  the  oscillating  part  whose 
time  average  is  zero.  (Following  Kaula  [1964],  we  write: 
k,  sin  elmpq,  though  notation  kfaimpq) sin  efXimpq)  would  be 
more  logical,  as  kj  and  e;  are  the  absolute  value  and  the  neg¬ 
ative  phase  of  the  complex-valued  function  kfximPq )■)  Here 
Flmp(i )  are  the  inclination  functions,  while  Gtpq(e )  are  the 
eccentricity  polynomials  identical  to  the  Hansen  coefficients 
Xg-~2 p+q)  2p) ■  The  Love  numbers  kj  are  our  |&/(x)|,  while 
the  phase  lags  are  given  by 

Qmpq  =  [(I  -  2p)u  +  (l-2p  +  q)M  +  m(Cl- e)\Atimpq, 

(36) 

where  9  is  the  sidereal  angle  of  the  primary,  9  is  its  spin  rate, 
A timpq  is  the  time  lag  at  the  mode 

uJbnpq  =  (l  -  2p)u)  +  {l-2p  +  q)M  +  m((l-8),  (37) 

while  a,  e,  i,  to,  f l,  M  stand  for  the  orbital  elements  of  the 
secondary:  the  semimajor  axis,  eccentricity,  inclination, 
argument  of  the  pericenter,  longitude  of  the  node,  and  mean 
anomaly.  The  actual,  physical  frequencies  of  deformation  in 
the  primary  are  given  by 

Ximpq  =  \uimpq\  =  \0  ~  2p)uj  +  (/  -  2p  +  q)M  +  m(Cl  -  e)\. 

(38) 

Below  we  shall  neglect  the  precessions  (il  ~  0,  u>  ~  0), 
and  shall  assume  that  M.  =  Ai0  +  n  ~  n. 

[50]  A  comprehensive  derivation  of  (35)— (38)  is  given 
by  Efroimsky  and  Williams  [2009].  As  each  term  in  the  sum 
(35)  is  proportional  to  ki  sin  e/mpq  =  \k^Ximpq)\sin  ei(xi„,pq)  = 
~Fm[k /(ximpq)]-  it  would  be  incorrect  tojntroduce  the  fre¬ 
quency  dependencies  of  sin  e/mpq  and  &/  =  \k/\  separately.  Both 


62200  =  922200  A?2200  —  X2200  At2200  Sgn  W2200-  (41) 

(We  ignore  the  mode  Impq  =  2000,  because  an  Impq  term  of 
the  torque  (35)  is  proportional  to  m.  Thence  the  said  mode 
contributes  nothing  in  the  spin  dynamics  of  the  tidally  dis¬ 
torted  body.  This  mode,  though,  does  influence  the  varia¬ 
tions  of  the  body’s  shape  and,  consequently,  the  dissipation 
rate.)  While  the  despinning  is  going  on,  i.e.,  while  9  >  n,  the 
tidal  mode  ca2 200  remains  negative,  as  can  be  seen  from  (39). 
Then  the  phase  lag  (41),  too,  remains  negative,  as  the  time 
lag  A?22oo  is  positively  defined.  Hence  the  imaginary  part  of 
the  Love  number,  T/w[A'2(X220o)]  =  ~k2  sin  £2200,  stays 
positive  during  the  despinning.  Naturally,  the  leading  term 
of  the  constant  part  of  the  torque, 

3 

T 2200  =  2  7  3Ts2ec  R5  a  6  k2  sin  £2200 
3 

=  ^  7  MsecR5a~6  k2  sin  |e220o|sgn  W2200,  (42) 

remains  negative  and  thereby  decelerating  the  spin.  Here  9  is 
the  primary’s  spin  rate,  while  Msec  is  the  mass  of  the  secondary. 
Recall  that  in  our  setting  the  role  of  the  tidally  distorted  primary 
is  played  by  the  satellite  (Iapetus),  while  the  role  of  the  tide¬ 
raising  secondary  is  played  by  the  planet  (Saturn). 

[52]  The  derivation  of  (42),  provided  by  Efroimsky  and 
Williams  [2009]  and  other  sources,  has  k2  sin  e22oo  approxi¬ 
mated  with  k2  Q22oo  sgn  tu22oo-  F°r  Iapetus,  we  cannot  afford 
this  approximation  because  over  some  period  of  its  history 
the  Q  factor  becomes  small  and  the  difference  between  Q  and 
sin|e|  becomes  considerable.  Thus,  we  have  to  stick  to  using 
^sin  £2200- 

[53]  Let  us  introduce  the  maximal  moment  of  inertia  of 
the  primary  (Iapetus): 

C  =  i  Mprim  R2 .  (43) 

Mprim  and  R  being  the  primary’s  mass  and  radius,  and  £  being 
a  numerical  factor  (equal  to  £  =  2/5  for  homogeneous 
spherical  bodies).  Then,  in  the  accepted  approximation,  the 
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despinning  rate  will  be  given  by 


3 

9  =  y  7  MXC  r3  a  6  sin  |e220o|sgn  UJ2200  +  0(e2sine) 

^  s  M-prim 

+  O (I2  sine) 


(44a 


or,  in  terms  of  the  mean  motion  n  =  yj 7  ( Mprim  +  Msec )  a~ 3 ,  by 


9  = 


n~Mi 


k2  sin  |e22oo|sgn  072200 
3  n2  Msex.  (R\ 3 


^  £>  Mprim  (Mprim  H“  Msecj  \Cl 

+  0(e2  sin  e)  +  0(f  sin  e)  =  ^  ^  ^ 

•  sin|e22oo|sgn  ut22oo  +  0(<?2sine)  +  0(;2sine) 

+  O  [Mprim / Msec) ,  (44b) 


where  we  recalled  that  Mprim  =  MIapetlls  <  Msec  =  MSatur„. 
As  follows  from  (39),  during  the  despinning  process  the  sign 
of  the  principal  tidal  mode  is  negative:  sgn  W2200  =  _1.  The 
value  of  sin|e22oo|  is  obtained  by  integration  over  layers,  as 
explained  in  section  3.5. 

[54]  In  the  approximation  of  a  homogeneous  primary  (and 
in  understanding  that  sin  e  =  e  +  0(e2)  and  sgn  u>  =  -1), 
expression  (44)  acquires  the  simple  form  of 

e=-(9  -  n)K.  +  0(e2e)  +  0(/2e)  +  0(e2).  (45) 


where 


£ 


7  M2 

sec 

C  Mpnm 


R2 

—T  k2  A?2200 

a° 


3 


n2M2 


C  Mprim  (Mprim  A  Msec ) 


3 

kl  A/2200  • 


(46) 


However,  this  simplification  is  deceptive  because,  to  imple¬ 
ment  it  in  practice,  one  still  has  to  calculate  the  frequency 
dependence  of  k2  At220o-  Based  on  the  dependence  of  the 
complex  Love  number  upon  the  complex  compliance,  such  a 
calculation  provides  the  answer  in  the  form  of  k2  sin  e22oo  = 
|^2(X220o)|sin  e2  (X2200)  =  _Tot[^2(X220o)]  anyway. 


5.  Some  Previous  Works  on  Iapetus’  Tidal 
Despinning 

5.1.  Peale  [1977]  and  Aleshkina  [2009] 

[55]  The  simplest  approach  to  despinning  is  to  assume  that 
both  k2  and  Q  are  constant.  Using  the  standardized  language 
of  the  empirical  power  law  Q  ~  (£  \f,  where  £  is  a  constant 
having  dimensions  of  time,  we  can  express  the  said  choice  as 

k2(x)  =  const,  (47  a) 


Q  ~  (£x)P  1  where  p  =  0.  (47b) 

[56]  This  way,  the  temperature  and  frequency  dependen¬ 
cies  and  a  lot  of  other  physics  get  ignored.  Within  this 
approach,  Peale  [1977]  explored  despinning  of  the  large 
moons  in  the  solar  system.  Setting  Q  =  100,  for  Iapetus,  Peale 
obtained  a  despinning  time  of  about  1 0  Ga,  which  was  much 
longer  than  for  other  large  satellites. 


[57]  Following  Peale  [1977],  Aleshkina  [2009]  set  the  Love 
number  constant.  In  her  dynamical  simulations,  Aleshkina 
[2009]  employed  the  so-called  modified  MacDonald  torque, 
i.e.,  the  MacDonald  [1964]  empirical  model  amended  with 
an  assumption  that  for  all  tidal  frequencies  x  the  appropriate 
time  lags  are  equal  to  the  same  quantity  £ : 

At  (X)  =  £,  (48) 

whence  the  tidal  Q  scales  as  inverse  frequency  \.  All  in  all: 

k2{x)  =  const,  (49a) 

Q  ~  (£x)Pi  where  p  =  —  1.  (49b) 

[58]  The  empirical  theory  of  bodily  tides  proposed  by 
Gerstenkorn  [1955]  and  furthered  by  MacDonald  [1964] 
suffered  an  inherent  inconsistency,  because  it  treated  the 
tidal  quality  factor  Q  as  a  constant.  The  empirical  theory 
becomes  equivalent  to  the  rigorous  development  by  Darwin 
[1879,  1880]  and  Kaula  [1964]  only  for  materials  whose 
time  lag  and  Q  factor  scale  as  (48)  and  (49b),  correspondingly 
[ Efroimsky  and  Williams,  2009].  From  this  viewpoint,  the 
use  of  (49b)  was  justified.  Unfortunately,  though,  (48)  and 
(49)  are  incompatible  with  the  actual  properties  of  solids.  This 
explains  why  employment  of  these  formulae  led  to  an 
implausibly  long  despinning. 

5.2.  Castillo -Rogez  et  al.  [2007] 

[59]  Castillo -Rogez  et  al.  [2007]  pointed  out  that  for 
Iapetus  to  despin  over  the  Solar  system  age,  its  quality  factor 
averaged  over  that  entire  time  span  must  have  been  less  than 
70  (for  k2  ~  10~3).  This  value  is  relatively  low.  To  obtain  it 
within  the  model  used  by  Castillo -Rogez  et  al.  [2007],  the 
internal  temperature  of  the  moon  ought  to  have  reached  at 
least  263  K,  i.e.,  to  have  approached  the  water-ice  melting 
point.  As  was  demonstrated  by  Castillo -Rogez  et  al.  [2007], 
this  might  be  achievable  due  to  the  insulating  effect  of 
porosity.  However,  a  high  porosity  is  compatible  neither  with 
freezing  of  the  fossil  equatorial  bulge,  nor  with  the  expected 
high  temperature  of  the  interior.  To  compensate  for  the  low 
porosity,  Castillo -Rogez  et  al.  [2007]  suggested  that  heating 
due  to  the  decay  of  short-lived  radioisotopes,  mainly  26A1, 
could  drive  early  compaction.  The  consequent  increase  in 
thermal  conductivity  would  then  justify  the  use  of  a  colder 
model  with  a  rapidly  thickening  lithosphere  and  with  the  deep 
interior  still  capable  of  reaching  temperatures  high  enough  to 
promote  dissipation  and  despinning. 

[60]  While  sufficiently  comprehensive,  the  combined 
geophysical  and  dynamical  calculation  by  Castillo  -Rogez 
et  al.  [2007]  suffered  two  drawbacks.  First,  it  relied  on  a 
simplistic  dissipation  model,  the  one  by  Maxwell.  Second, 
Castillo -Rogez  et  al.  [2007]  did  not  include  convection.  As 
convection  is  initiated,  it  prevents  further  warming  of  the 
interior.  This  problem  was  tackled  later  on  by  Robuchon  et  al. 
[2010], 

5.3.  Robuchon  et  al.  [2010] 

[61]  Keeping  some  of  the  assumptions  used  by  Castillo- 
Rogez  et  al.  [2007]  and  introducing  two  major  upgrades, 
Robuchon  et  al.  [2010]  recently  developed  a  new  model  of 
Iapetus.  The  first  upgrade  was  a  block  modeling  convective 
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Maxwell  Body  Kelvin-Voigt  Element 

It  r'®®®5'— 1 


Figure  1.  Schematic  representation  of  various  viscoelastic 
and  anelastic  models.  A  model  is  illustrated  by  an  assemblage 
of  springs  and  dashpots.  Springs  (labeled  with  p )  are  aimed  to 
represent  the  elastic  properties  of  the  medium,  while  dashpots 
(denoted  with  ?y)  stand  for  the  viscous  properties.  In  the  Bur¬ 
gers  model,  anelasticy  is  represented  by  one  Voigt  element, 
which  is  characterized  by  its  specific  elastic  modulus  and  vis¬ 
cosity.  In  the  case  of  the  Andrade  model,  the  anelastic  com¬ 
ponent  represents  an  infinite  number  of  dashpots  in  series 
in  parallel  with  an  infinite  number  of  springs  and  represents 
a  continuous  distribution  of  compliances  and  thus  relaxation 
times. 


heat  transfer.  They  demonstrated  that  convection  can  start  at 
temperatures  lower  than  the  optimum  temperature  necessary 
to  achieve  significant  dissipation  in  the  Maxwell  model.  As  a 
second  upgrade,  they  resorted  to  a  different  dissipation 
model,  the  Burgers  element,  in  order  to  determine  whether  or 
not  a  dissipation  rate  leading  to  despinning  can  be  achieved 
at  low  temperature.  The  Burgers  element,  illustrated  by 
Figure  1,  should  not  to  be  confused  with  the  so-called 
extended  Burgers  model,  which  presents  a  more  complex 
distribution  of  relaxation  times.  Resorting  to  the  Burgers 
element  was  based  on  the  observation  by  Reeh  et  al.  [2003] 
that  the  model  can  reproduce  the  bending  of  floating  gla¬ 
ciers  in  response  to  ocean  tides.  Applying  the  model  to 
lapetus,  Robuchon  et  al.  [2010]  tested  different  values  for  the 
viscous  and  elastic  moduli  characterizing  the  transient 
response  (between  elastic  behavior  and  viscous  creep)  in 
order  to  find  a  combination  that  could  lead  to  a  dissipative 
interior  despite  convection  onset.  They  did  find  successful 
models  in  which  dissipation  could  occur  at  temperatures  as 
low  as  200  K  and  in  few  hundred  million  years.  This  way, 
Robuchon  et  aids  [2010]  study  was  a  pioneer  attempt  to  take 
into  account  the  anelastic  response  of  ice  in  a  tidal  dissipation 
model.  At  the  same  time,  it  should  be  pointed  out  that  the 
choice  of  the  Burgers  model  by  Robuchon  et  al.  [2010]  suf¬ 
fers  the  same  flaws  as  the  choice  of  the  Maxwell  model. 
Indeed,  there  exists  a  considerable  discrepancy  between  the 
Burgers  element  and  experimental  results  [Jellinek  and  Brill, 


1956;  Tatibouet  et  al.,  1987;  Jackson,  1993;  Cole,  1995; 
Cooper,  2002] .  The  reason  for  this  is  that  attenuation  is  driven 
by  the  motion  of  defects  whose  geometry  and  distribution 
are  complex  and  cannot  be  accounted  for  by  a  single  Voigt 
element. 

[62]  We  do  not  question  the  fact  that  the  Burgers  model  did 
provide  a  good  match  to  the  observations  of  Reeh  et  al. 
[2003],  However,  as  pointed  out  by  Sold  and  Spohn  [1997], 
who  aimed  to  interpret  Mars’  tidal  Love  number  and  dissi¬ 
pation  factor  at  the  period  of  the  forcing  frequency  exerted  by 
Phobos,  both  the  Maxwell  model  and  the  Burgers  element  can 
yield  a  result  consistent  with  that  observation.  However,  these 
different  models  imply  different  values  for  the  viscosity  of 
Mars’  mantle.  Thus,  in  principle,  a  single  observational 
constraint  may  be  interpreted  via  several  possible  dissipation 
models,  which  questions  the  geophysical  meaning  of  the 
input  parameters  to  these  models. 

[63]  Besides,  forward  modeling  of  planetary  dissipation 
involves  a  different  set  of  issues.  For  one,  the  input  param¬ 
eters  to  the  Burgers  model  (for  the  transient  creep)  inferred  for 
Earth’s  cryosphere  or  mantle  are  not  necessarily  applicable  to 
other  bodies.  The  values  inferred  for  Earth  are  for  specific 
conditions  of  stress  and  temperature  that  are  far  different  from 
those  expected  in  lapetus  or  other  icy  satellites:  a  temperature 
of  -30°C  in  the  case  of  the  Reeh  et  al.  [2003]  study  [after 
Brill  and  Camp,  1961]  and  also  a  stress  of  a  few  hundred 
kilopascals.  Thus,  the  extrapolation  of  the  Burgers  model  to 
planetary  bodies  is  not  straightforward.  Robuchon  et  al. 
[2010]  varied  the  possible  parameters  of  the  transient 
Burgers  element  before  finding  the  combinations  that  lead  to 
despinning.  As  explained  in  section  6,  the  transient  and  steady 
state  responses  are  necessarily  coupled  since,  in  the  low- 
stress  conditions  relevant  to  lapetus’  despinning,  they  both 
involve  the  same  defects. 

[64]  As  a  summary,  since  the  experimental  data  on  atten¬ 
uation  in  silicates  and  other  materials  do  not  fit  the  Burgers 
element  well,  we  would  use  caution  when  applying  that 
model. 

6.  What  Tidal  Friction  Model  for  lapetus? 

6.1.  On  the  Limitations  of  “Tin  Toys  Models” 

[65]  Numerous  studies  assume  planetary  materials  to 
behave  as  a  Maxwell  body,  a  model  represented  by  a  purely 
viscous  damper  and  a  purely  elastic  spring  connected  in 
series.  The  model  was  applied,  for  example,  to  a  variety  of  icy 
satellites  [e.g.,  Tobieetal.,  2005;  Castillo -Rogezetal.,  2007], 
The  Maxwell  model  entails  the  following  expression  for  the 
phase  lag: 


tan  6{x) 


1 

X 


(50) 


X  being  the  frequency.  Being  unable  to  account  for  the 
anelastic  behavior  observed  in  planetary  materials  [Cooper, 
2002],  the  model  fails  to  account  correctly  for  attenuation 
in  minerals  (see  Karato  [2008]  for  a  detailed  study,  or 
Efroimsky  and  Lainev  [2007]  for  a  very  brief  review). 

[66]  As  discussed  below  in  more  detail,  numerous  labora¬ 
tory  measurements  indicate  that  the  transient  response  of 
planetary  materials  is  complex,  and  cannot  be  simply 
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accounted  for  by  a  discrete  set  of  relaxation  peaks.  An  ade¬ 
quate  description  requires  a  complicated  distribution  function 
of  the  relaxation  time  [e.g.,  Jackson  et  al.,  2002;  Cooper, 
2002],  The  reason  for  broadening  of  the  relaxation  peak, 
compared  to  the  Debye  peaks  characterizing  the  Maxwell 
body  or  the  Burgers  element,  is  that  in  realistic  materials  the 
defects  driving  dissipation  (e.g.,  grain  boundaries  or  dis¬ 
locations)  are  not  homogeneously  distributed  in  the  material. 
Besides  this,  the  defects’  parameters  (e.g.,  dislocations 
lengths  or  grain  shapes)  may  vary  by  orders  of  magnitude 
within  one  icy  shell.  These  circumstances  complicate  the 
picture,  and  make  simplistic  tin  toys  models  (i.e.,  assem¬ 
blages  of  springs  and  dashpots)  insufficient. 

6.2.  The  “Elbow”  Dependence  and  the  Andrade  Model 

[67]  Experimental  exploration  of  attenuation  in  planetary 
materials  began  in  the  middle  of  the  20th  century.  In  appli¬ 
cation  to  minerals,  this  research  was  motivated  by  the  need  to 
interpret  seismic  data.  In  application  to  ices,  it  was  motivated 
by  the  interest  in  the  mechanical  properties  of  ice  shelves;  so 
the  pioneering  data  were  obtained  from  observing  the  tran¬ 
sition  from  elastic  deformation  to  steady  state  creep  (e.  g.,  to 
the  primary  creep  of  ice  under  a  constant  load).  We  refer  the 
reader  to  McCarthy  and  Castillo  -Rogez  [201 1]  for  details  on 
the  history  of  research  on  the  dissipative  properties  in  ices. 

[68]  These  measurements  demonstrate  that  the  slope  of 
the  attenuation  spectrum  undergoes  a  drastic  change  at  a 
threshold  frequency  that  marks  the  boundary  between  a  low- 
frequency  and  high-frequency  regime  (see  Efroimsky  and 
Lainey  [2007]  for  a  review). 

[69]  This  entails,  as  a  possible  option,  the  following  fre¬ 
quency  dependence: 

QJl  =  tan  S  =  (£  x)p,  where  p  =  0.2  —  0.4 

for  X»Xo  and  p  =  ~  1.0  for  x  <  Xo,  (51) 

The  empirical  parameter  £  has  the  dimensions  of  time,  and 
has  the  physical  meaning  of  the  average  timescale  associated 
with  the  dominating  dissipation  mechanism.  As  different 
mechanisms  may  dominate  over  different  ranges  of  fre¬ 
quencies,  both  £  and  p  can  bear  dependence  on  %,  but  these 
dependencies  are  assumed  to  be  slow,  except  for  the  fast 
change  of  p  at  the  threshold  frequency  \  =  Xo- 

[70]  Evidence  for  the  applicability  of  (51)  at  high  fre¬ 
quencies  comes  from  measurements  in  a  lab  [e.g.,  Gribb  and 
Cooper,  1998;  Jackson  et  al.,  2002]  and  seismic  data  [Shito 
et  al.,  2004],  For  ices,  such  behavior  has  been  reported,  for 
example,  by  Cole  [1990].  Observational  evidence  for  the 
low-frequency  band  comes  from  observation  of  the  Chandler 
wobble  and  of  the  mantle’s  response  to  postglacial  rebound 
[Romanowicz  and  Durek,  2000;  Karato  and  Spetzler,  1990], 
The  change  in  the  frequency  dependence  indicates  that  the 
attenuation  in  the  material  is  dominated  at  low  frequencies  by 
its  viscous  properties  [Karato,  2008],  while  at  higher  fre¬ 
quencies  it  is  dominated  by  anelasticity  [e.g.,  Jackson  et  al., 
2002], 

[71]  Several  experimental  studies  have  demonstrated  that 
the  change  in  regime  is  not  immediate  but  involves  a  transi¬ 
tional  region  [e.g.,  Jackson  et  al. ,  2002],  Jackson  et  al.  [2002] 
also  demonstrated  on  theoretical  grounds  that  in  the  case  of 


grain  boundary  diffusion,  the  adjustment  between  the  two 
regimes  occurs  on  a  timescale  that  is  a  function  of  the 
microstructural  and  viscoelastic  properties  of  the  material. 
However,  further  experimental  and  theoretical  studies  are 
needed  in  order  to  characterize  the  threshold  xo- 

[72]  Analysis  of  numerous  forced  oscillation  experiments 
indicate  that  the  models,  which  include  a  broad  relaxation 
term  (such  as  the  Andrade  model  or  the  extended  Burgers 
model),  provide  the  best  fit  to  experimental  data  [e.g.,  Tan 
et  al.,  2001;  Jackson  et  al.,  2002;  Cooper,  2002]. 

[73]  Over  the  past  years,  the  model  pioneered  by  Andrade 
[1910,  1914]  has  been  reported  match  for  a  variety  of  mate¬ 
rials,  over  a  wide  range  of  experimental  conditions  [Cottrell 
and  Aytakin,  1947;  Duval,  1978].  The  abundance  of  experi¬ 
mental  data  supporting  the  model  is  the  reason  why  we  made 
it  our  choice.  The  model  adequately  describes  the  nonlinear, 
transient  response  both  in  metals  [Andrade,  1914]  and 
minerals  [e.g.,  Jackson,  1993]. 

[74]  In  the  time  domain,  the  creep  function  corresponding 
to  the  Andrade  model  is 

J(t-t')  =  JQ(t-t')  +  f3(t  —  t')a  +  -  (t  —  t'),  (52) 

v 

r]  denoting  the  steady  state  viscosity,  and  J=  J(0)  =  1/^(0)  = 
1/p,  being  the  unrelaxed  compliance,  while  a  and  (3  standing 
for  empirical  parameters.  Here  we  employ  the  Heaviside 
function  0(7  -  t')  to  make  sure  that  differentiation  of  J(t  -  t'), 
with  subsequent  multiplication  by  ffyit')  and  integration  over 
t',  generates  the  right  expression  for  the  strain,  with  the 
unrelaxed  term  Jo^v  present. 

[75]  Parameter  0  characterizes  the  intensity  of  anelastic 
friction  in  the  material,  and  therefore  must  depend  upon  the 
density  of  the  defects.  The  shape  of  this  dependence  remains 
unknown,  because  no  research  has  ever  been  undertaken  in 
this  direction  in  the  case  of  diffusion-driven  attenuation. 
(In  this  regard,  it  would  be  interesting  to  mention  the 
anelastic -viscoelastic  model  developed  by  Cole  [1995],  who 
managed  to  express  the  complex  compliance  as  an  explicit 
function  of  the  defect  density  in  the  case  of  dislocation-driven 
attenuation.)  Experimental  measurements  [e.g.,  Tan  et  al., 
2001;  Jackson  et  al.,  2002]  show  that  the  value  of  0 
depends  upon  the  temperature.  Theoretical  work  summarized 
by  Karato  and  Spetzler  [1990]  also  suggests  that  the  attenu¬ 
ation  amplitude  should  be  a  function  of  chemistry  (e.g., 
fugacities).  Parameter  a  determines  the  duration  of  the  tran¬ 
sient  response  in  the  primary  creep.  It  depends  upon  the  stress 
and  upon  the  relaxation  time  of  the  dominating  mechanism  of 
friction  [Castelnau  et  al.,  2008].  The  values  of  a  for  olivine- 
rich  rocks  typically  fall  within  the  interval  0. 1-0.5,  most  often 
within  0. 2-0.4.  A  remarkable  experimental  fact  is  that  the 
water  ice,  despite  all  its  physical  and  chemical  differences 
from  minerals,  obeys  this  same  law,  with  the  parameter  a 
having  values  similar  to  those  it  has  for  rocks.  This  has  been 
indicated  by  experiments  carried  out  in  the  grain  boundary 
sliding  and  in  the  dislocation  creep  regimes  [Glen,  1955; 
Jellinek  and  Brill,  1956;  Duval,  1978;  Cole,  1995;  Castelnau 
et  al.,  2008]. 

[76]  The  Andrade  model  can  be  illustrated  with  a  Maxwell 
element  standing  in  series  with  an  infinite  series  of  dashpots, 
set  in  parallel  with  an  infinite  series  of  springs  (Figure  1 ).  This 
configuration  reproduces  the  distribution  of  relaxation  times 
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Figure  2.  Dissipation  spectrum  for  the  Andrade  model  (thick  line)  computed  from  equation  (54),  with 
input  parameters  assuming  values  relevant  to  the  case  of  Iapetus:  a  =  0.33,  ft  =  10~12,  p  =  3.3  GPa,  and  r]  = 
1015  Pa  s.  The  dashed  curve  represents  the  corresponding  Maxwell  model  computed  for  the  same  values  of 
unrelaxed  shear  modulus  and  viscosity.  The  thick  grey  curve  represents  the  slope  of  the  spectrum  computed 
with  the  Andrade  model.  The  boxes  with  arrows  indicate  the  trend  followed  by  this  slope  for  infinitely  high 
or  infinitely  low  frequencies. 


characterizing  the  complex  transient  response  of  the  medium. 
This  peak-broadening  construction  is  inherent  also  in  some 
other  models  similar  to  that  of  Andrade  [e.g.,  Cole,  1995], 
As  all  the  other  empirical  models,  the  Andrade  model  is 
“parametrically  economical,”  though  it  also  shares  the  main 
drawback  of  such  models,  the  “lack  of  physical  transparency” 
[Jackson,  1993].  However,  in  recent  years  the  transient  creep 
described  by  the  Andrade  model  was  successfully  derived 
via  more  fundamental  approaches  describing  the  motion  of 
defects  under  stress  for  different  types  of  defects,  grain 
boundaries  in  ice  [Castelnau  et  al.,  2008],  or  dislocation 
jamming  in  ice  [. Miguel  et  al.,  2000]. 

[77]  The  complex  compliance  is  inferred  from  the  creep 
function  (52)  through  the  Laplace  or  Fourier  transform,  pro¬ 
vided  that  the  material  behaves  linearly  (which  is  surely  true 
under  weak  forcing).  According  to  Findley  et  al.  [1976]  and 
Gribb  and  Cooper  [1998]  the  complex  compliance  for  a 
material  obeying  the  Andrade  model  is  given  by 

J(x)  =  J-  —  +  mra  r  (l+a),  (53) 

vx 

J  =  Up  being  the  unrelaxed  compliance,  \  denoting  the  fre¬ 
quency,  and  T  standing  for  the  Gamma  function.  Dependence 


(53)  entails  the  following  expression  for  the  phase  lag  as  a 
function  of  the  frequency  [e.g.,  Nimmo,  2008]: 

t-_  c  _  ImQ)  _  (hxr1  +  X~a  ft  sin(f)r(a  +  1) 

Re(J)  J  +  jc*  p  cos(^)  r(a  +  1)  ’  1  j 

whence  we  see  that 


tan <5  ~  x  Q  forx>Xo, 

tan<5~x^  for  x  <  Xo- 


(55) 


[78]  In  section  6.3.3,  we  explain  that  the  parameter  ft  can 
be  expressed  in  a  simple  way  (63)  via  J  and  p.  Using  that 
interrelation,  it  will  be  easy  to  demonstrate  that  the  threshold 
frequency  is  the  inverse  Maxwell  time: 

Xo  =  tm'  =  ^  =  (J  r])~\  (56) 


[79]  The  attenuation  spectrum  derived  theoretically  from 
equation  (54),  with  the  parameters’  values  realistic  for  our 
study,  is  presented  in  Figure  2.  Besides  the  disagreement 
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between  the  Andrade  and  Maxwell  model-based  curves  at 
high  frequencies,  which  is  well  documented  by  laboratory 
experiments,  Figure  2  shows  also  a  difference  between  the 
models  at  very  low  frequencies.  With  the  insufficient  obser¬ 
vational  data  available  at  low  frequencies,  it  is  impossible  to 
say  which  of  the  two  curves  is  a  better  fit  over  that  band.  Some 
scholars  favor  the  Maxwell  model,  assuming  that  the  mantle 
is  strictly  viscous  in  the  low-frequency  limit  (S.  Karato, 
personal  communication,  2010).  We  shall  approach  this  sit¬ 
uation  in  section  6.4  by  bringing  in  a  composite  scaling  law, 
which  employs  the  Andrade  model  to  describe  the  anelas- 
ticity  of  the  material  in  the  high-frequency  regime  and 
switches  to  the  Maxwell  model  at  low  frequencies. 

6.3.  Mechanisms  Driving  Internal  Friction  in  Iapetus’ 
Conditions 

[so]  We  open  this  section  with  presenting  the  mechanisms 
that  are  likely  to  drive  attenuation  in  Iapetus  (sections  6.3.1- 
6.3.3).  Then  we  describe  how  the  corresponding  physics 
was  implemented  in  the  geophysical  model  (section  6.3.4). 
Finally,  we  illustrate  these  principles  with  several  numerical 
runs  (section  6.3.5). 

6.3.1.  Dominating  Mechanisms 

[si]  Ice  creep  involves  several  mechanisms  whose  activa¬ 
tions  E*  energies  are  close  (50-60  kJ  moF1  [ Goldsbv  and 
Kohlstedt,  2001]).  The  creep  rate  is  expressed  as  a  function 
of  the  contribution  of  each  mechanism  [e.g.,  Goldsby  and 
Kohlstedt ,  2001]: 

u  =  ilBD  +  Mvd  +  «dc  +  (—  +T  J  (57) 

\IIGBS  UBS  J 

where  BD  refers  to  the  boundary  diffusion  creep  (also  called 
Coble  creep),  VD  refers  to  the  grain  volume  diffusion  creep 
(also  called  self-lattice  diffusion  orNabarro-Herring  creep,  or 
self-diffusion,  or  bulk  diffusion),  DC  refers  to  the  dislocation 
creep,  GBS  refers  to  the  grain  boundary  sliding,  and  BS  refers 
to  the  basal  slip.  The  relative  contributions  from  these 
mechanisms  depend  on  the  material  grain  size  (for  BD,  VD, 
and  GBS)  and  the  magnitude  of  the  stress  (for  DC,  BS 
and  GBS).  The  stress  magnitude  determines  the  amount  of 
mechanical  energy  available  to  activate  the  motion  of  defects 
or,  in  the  case  of  diffusion,  to  activate  breaking  of  bonds  and 
reorientation  of  protons. 

[82]  In  Iapetus,  the  tidal  stress  is  expected  to  be  weak. 
Being  of  the  order  of  1 0  4  MPa,  it  is  at  least  three  orders  of 
magnitude  lower  than  the  tidal  stress  undergone  by  Encela- 
dus.  Accordingly,  tidal  displacements  of  Iapetus’  surface  are 
about  2  cm,  while  the  strain  rate  does  not  exceed  1CT13  s-1. 
These  values  of  the  tidal  displacement  and  strain  are  of  the 
same  orders  of  magnitude  as  the  values  of  the  displacement 
and  strain  generated  by  the  convective  stress  [e.g.,  Durham 
and  Stern,  2001].  Accordingly,  for  low  tidal  stressing,  one 
should  expect  dissipation  in  Iapetus  to  be  dominated  by  the 
diffusion  creep  [ Goodman  et  al.,  1981;  Goldsby  and 
Kohlstedt,  2001;  Barr  and  Pappalardo,  2005]. 

[83]  Diffusion  creep  is  known  to  involve  motion  of  both 
ionic  and  orientational  (Bjerrum)  defects  [Onsager  and 
Runnels,  1969].  The  latter  type  of  defects  has  been  identi¬ 
fied  as  an  important  source  of  mechanical  and  electrical 
relaxation.  Ice  is  one  of  the  very  few  materials  in  which  not 


only  translational  motion  but  also  rotational  motion  of  defects 
takes  place  under  stress.  Translational  is  the  self-lattice  dif¬ 
fusion,  i.e.,  diffusion  of  material  from  regions  under  loading 
to  regions  under  tension.  Rotational  is  reorientation  of  the 
Bjerrum  defects  (or  spin  lattice  diffusion),  which  has  been 
much  documented  due  to  its  notable  signature  in  ice  dielectric 
properties.  Another  mechanism,  the  Coble  [1963]  creep,  is 
diffusion  of  interstitials  at  the  grain  boundaries.  Being  largely 
responsible  for  seismic  wave  attenuation  in  the  Earth’s  mantle 
[e.g.,  Gribb  and  Cooper,  1998],  this  process  is  of  a  lesser 
importance  in  cold  pure  water  ice  (see  below  for  details).  For 
the  stress  magnitude  and  the  grain  size  considered  in  this 
study,  diffusion  may  be  limited  by  the  basal  slip  at  tempera¬ 
tures  below  120  K  [Goldsby  and  Kohlstedt,  2001].  At  these 
low  temperatures,  defect  motion  is  limited  (the  strain  rate  is 
lower  than  10~29  s  *),  so  that  the  contribution  of  basal  slip  to 
the  overall  dissipation  history  of  Iapetus  is  negligible. 

[84]  In  the  following  sections  we  study  the  contributions  of 
various  types  of  defect  motion  in  the  ice  lattice,  assuming  that 
Iapetus  is  made  up  of  a  pure  water  ice,  an  approximation 
conventionally  used  for  satellites  that  contain  a  small  fraction 
of  silicates. 

6.3.2.  Proton  Reorientation 

[85]  Stress-induced  reorientation  of  water  molecules, 
assisted  by  the  diffusion  of  orientational  (Bjerrum)  defects, 
results  in  a  small  anelastic  strain  whose  characteristic  relax¬ 
ation  time  is  strongly  dependent  on  temperature  and  fre¬ 
quency  and  is  very  close  to  the  dielectric  relaxation  time 
[Petrenko  and  Whitworth,  1999].  The  temperature  range, 
over  which  proton  reorientation  is  most  effective,  is  defined 
by  the  forcing  frequency.  For  the  frequencies  experienced  by 
Iapetus  over  its  history,  this  range  is  from  about  150  K  down 
to  125  K.  Therefore  we  have  to  include  proton  reorientation 
in  our  study,  as  a  mechanism  dominating  friction  at  these 
temperatures. 

[86]  Proton  reorientation  is  independent  of  the  grain  size  or 
stress  [Tatibouet  et  al.,  1983],  Although  the  phase  lag  mea¬ 
surements  by  Tatibouet  et  al.  [1981,  1983]  provide  the  only 
available  data  on  proton  reorientation  at  excitation  frequen¬ 
cies  below  1(T3  Hz,  Petrenko  and  Whitworth  [1999]  have 
confirmed  the  consistency  of  these  results  with  electric 
measurements.  Besides,  the  increase  in  dissipation  measured 
at  temperatures  from  167  K  to  130  K  is  consistent  with  the 
increased  mobility  of  protons  [Johnson  and  Quickenden, 
1997]. 

[87]  Mechanical  measurements  indicate  that  the  relaxation 
peak  temperature  (the  temperature  at  which  the  dissipation 
rate  is  maximal)  decreases  with  decreasing  frequency:  from 
167  K  at  1  Hz  to  150  K  at  1 0  Hz,  as  reported  by  Tatibouet 
et  al.  [1981].  For  frequency  10  4  Hz  and  temperature  130  K, 
these  authors  obtained  tan<5  =  0.02.  Extrapolating  the  results 
linearly  to  lower  temperatures,  they  suggested  that  the  quality 
factor  could  be  as  low  as  40  at  the  temperature  of  125  K  and 
the  frequency  of  4  x  1 0  3  Hz  (corresponding  to  tidal  excita¬ 
tion  of  Iapetus  with  a  spin  period  of  10  h).  However,  their 
extrapolation  is  not  consistent  with  the  fact  that  the  proton 
mobility  starts  decreasing  beginning  from  125  K  and 
becomes  very  limited  at  100  K  [Johnson  and  Quickenden, 
1997],  which  most  probably  entails  a  decrease  in  dissipation 
rate.  In  any  case,  these  measurements  indicate  that  between 
100  and  167  K,  the  quality  factor  associated  with  proton 
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Figure  3.  Attenuation  resulting  from  proton  reorientation  in  ice,  as  a  function  of  the  temperature, 
computed  using  the  parameters  from  Tatibouet  et  al.  [1981],  The  temperature  of  the  relaxation  peaks 
is  shown  as  a  function  of  frequency.  The  calculation  is  carried  out  for  the  temperatures  down  to  120  K, 
below  which  proton  mobility  gets  significantly  decreased.  For  details,  see  section  6.3.2,  specifically 
formulae  (58)-(60). 


reorientation  is  lower  than  or  close  to  1 00.  The  relaxation  time 
for  that  mechanism  is  given  by  [Vassoille  et  al,  1974; 
Tatibouet  et  al.,  1981]: 


sxp(Em/RT) 
c0  +  exp  (-Ef/RT) 


(58) 


where  Ef(~. 25  kJ  moF1)  and  Em{~ 35  kJ  mol-1)  are  the 
Bjerrum  defects’  energies  of  formation  and  migration, 
respectively;  while  Co  is  the  concentration  of  the  extrinsic 
defects  that  determine  the  formation  of  the  Bjerrum  defects. 
As  suggested  by  Tatibouet  et  al.  [1981],  we  take  r0  ~  6  x 
1CF16  s.  Tatibouet  et  al.  [1981]  also  demonstrated  that  the 
value  of  Co  lies  between  5  x  1CF8  and  1CF7.  From  their 
experimental  measurements,  they  established  an  empirical 
formula  for  the  maximal  phase  lag  (the  one  at  the  dissipation 
peak)  as  a  function  of  temperature: 

tan<Smax=||,  (59) 

with  Ep  ~  25  kJ  mol  1 . 

[88]  Using  the  available  phase  lag  measurements  and  tak¬ 
ing  into  account  that  the  friction  coefficient  has  a  Debye  peak 
(as  has  been  observed  by  electric  experiments),  one  can 
hypothesize  that  the  dissipative  behavior  dominated  by  pro¬ 
ton  reorientation  should  fit  the  S  AS  model  because  this  model 
permits  for  a  Debye  peak.  The  SAS  model  connects  the  stress 
and  the  strain  via  <7f„  +  T^oy,,  =  2/z(wf„  +  r„wtv)  with 


Ta  and  tu  being  characteristic  times,  and  dot  standing  for  a 
time  derivative.  In  the  frequency  domain,  this  entails:  2/2  = 
v  =  2q  |  +  wherefrom  we  obtain:  ta nb  =  Tm[jl]/ 

7 Ze[p]  =  [+~T°  ji-  Thence  it  is  easy  to  show  that  the  tangent  is 
related  to  its  maximal  value  through 

tan  8  =  2  [tan<5]max  ^  +T  , 

where  the  effective  relaxation  time  is  defined  as  f  =  yjTaTu. 
It  can  be  identified  with  the  time  tp  given  by  (58)  and  can  be 
experimentally  determined  using  the  observed  shape  of  the 
Debye  peak.  Combining  the  above  equations,  we  obtain  tan<5 
as  a  function  of  the  temperature  and  frequency  (see  Figure  3). 
Figure  3  shows  that  proton  reorientation  is  expected  to  be  an 
important  source  of  dissipation  at  high  forcing  frequencies 
and  low  temperatures,  i.e.,  at  the  conditions  whereto  Iapetus 
is  subject  following  accretion.  As  calculated  in  detail  below, 
this  mechanism  turns  out  to  be  the  dominant  source  of  dis¬ 
sipation  over  the  first  few  hundred  million  years  following  the 
formation  of  the  satellite. 

[89]  Finally,  note  that  we  lack  information  on  the  effect  of 

proton  reorientation  upon  the  effective  elastic  modulus.  So 
we  assume  that  the  latter  remains  constant  and  is  equal  to  the 
unrelaxed  shear  modulus  /i  =  0). 

6.3.3.  Diffusion-Assisted  Grain  Boundary  Sliding 

[90]  Although  it  is  established  that  diffusion  is  the  domi¬ 
nant  mechanism  responsible  for  the  Newtonian  viscosity  of 
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Table  1.  Input  Parameters 


Parameter 

Definition 

Value 

Source 

a 

Semimajor  axis  of  the  secondary 

3,560,820  km 

http://ssd.jpl.nasa.gov 

A,  B,  C 

Moments  of  inertia  of  the  satellite  (A  <  B  <  C) 

function  of  evolution 

Co 

Initial  concentration  in  Bjerrum  defects 

10-7 

Tatibouet  et  al.  [1981] 

d 

Ice  grain  size 

0.1  to  1  mm 

E>b,0 

Preexponential  term  for  the  grain  boundary  diffusion  coefficient 

6  x  nr4  m2  s-1 

Goldsby  and  Kohlstedt  [2001] 

A,0 

Preexponential  term  for  the  volume  diffusion  coefficient 

9.1  x  l(T4  nr  s-1 

Ramseier  [1967] 

e 

Eccentricity  of  the  secondary 

0.028  612  5 

http://ssd.jpl.nasa.gov 

E„ 

Activation  energy  for  grain  boundary  diffusion  creep 

49  kJ  moF1 

Goldsby  and  Kohlstedt  [2001] 

Ed 

Activation  energy  for  volume  diffusion  creep 

59.4  kJ  moF1 

Goldsby  and  Kohlstedt  [2001] 

Ef 

Activation  energy  of  Bjerrum  defect  formation 

25  kJ  moF1 

Tatibouet  et  al.  [1981] 

Em 

Activation  energy  of  Bjerrum  defect  migration 

35  kJ  mol  1 

Tatibouet  et  al.  [1981] 

Ep 

Activation  energy  for  proton  reorientation 

25  kJ  mol  1 

Tatibouet  et  al.  [1981] 

g 

Surface  gravity 

0.223  m  s  2  at  the  equator 

Msec 

Mass  of  the  secondary 

1.805  635  x  1021  kg 

Jacobson  et  al.  [2006] 

n 

Mean  motion 

R 

Satellite’s  mean  radius 

735.6  km 

Thomas  et  al.  [2007] 

Tm 

Melting  temperature 

273  K  (268  K  if  depressed  by  impurities) 

v„, 

Molar  volume  of  defects 

1.97  x  10~5  m3 

Fletcher  [1970] 

w 

Grain  boundary  thickness 

9.04  x  10~10  m 

Goldsby  and  Kohlstedt  [2001] 

a 

Input  parameters  to  the  Andrade  creep  equation 

0.2  to  0.5 

7 

Newton’s  gravitational  constant 

6.674  X  10~"  m3  kg~'  s~2 

R,  M(0) 

Unrelaxed  shear  modulus 

3.3  GPa 

Parameswaran  [1987] 

To 

Reference  relaxation  time  for  proton  reorientation 

6  x  10~16  s 

Tatibouet  et  al.  [1981] 

ice  at  low  stress  [e.g.,  Bromer  and  Kingeiy,  1 968],  it  has  been 
little  studied  experimentally,  because  this  phenomenon  is  not 
easily  achieved  in  laboratory  conditions  [Duval,  1978], 
Especially,  the  transient  response  of  ice  in  the  regime  of  self- 
diffusion  is  poorly  constrained.  The  following  description  is 
very  much  inspired  by  the  literature  on  silicates,  with  the 
obvious  caveat  that  experimental  measurements  are  neces¬ 
sary  to  confirm  this  approach  in  application  to  ices. 

[91]  The  strain  rate  corresponding  to  diffusion  creep  is 
[after  Goodman  et  al.,  1981;  Ban •  and  Pappalardo,  2005]: 


it 


42  Vm(7(ij 

iRoTind2 


(60) 


where  Vm  denotes  the  molar  volume  of  defects,  R  is  the  per¬ 
fect  gas  constant,  Tm  is  the  melting  temperature,  d  stands  for 
the  grain  size,  w  is  the  grain  boundary  thickness,  while  Dv  and 
Dh  are  the  diffusivities  of  self-diffusion  and  grain  boundary 
diffusion,  respectively.  According  to  Goldsbv  and  Kohlstedt 
[2001],  both  diffusivities  scale  with  temperature  as 

Dk  =  Dkfi  exp  (61) 


the  subscript  k  =  v,  b  referring  to  the  volume  diffusion  and 
grain  boundary  diffusion  mechanisms,  respectively.  The 
parameters  Dk  0  and  Ek  represent  the  preexponential  magni¬ 
tude  and  activation  energy  of  these  mechanisms.  The  values 
of  these  parameters  are  given  in  Table  1. 

[92]  In  a  pure  water  ice  at  temperatures  well  below  the 
premelting  regime,  the  term  ~\'Dh  is  several  orders  of  mag¬ 
nitude  smaller  than  Dv,  so  grain  boundary  diffusion  con¬ 
tributes  little  to  deformation  of  the  material.  However,  grain 
boundary  diffusion  may  become  significant  in  the  presence  of 
partial  melt  [Goldsbv  and  Kohlstedt,  2001]  or  impurities 
[Bromer  and  Kingerv,  1968],  We  explore  these  scenarios  in 
section  6.4. 


[93]  As  the  temperature  dependence  of  the  diffusivity  fol¬ 
lows  the  Arrhenius’  law,  it  is  possible  to  rewrite  (60)  as 

A  /— £v*\ 

"  =  ^expUr>  (62) 

u  signifying  the  shear  strain,  E*  standing  for  the  activation 
energy. 

[94]  Robuchon  et  al.  [2010]  have  suggested  that  diffusion- 
creep-driven  convection  in  lapetus  should  involve  ice  grains 
smaller  than  100  pm,  as  the  presence  of  impurities  prevents 
grain  growth.  During  compaction,  grains  undergo  commi¬ 
nution  (size  reduction)  due  to  brittle  fracture  [e.g.,  Durham 
et  al.,  2005].  So  we  take  the  said  value  as  an  upper  bound. 
We  set  the  lower  bound  at  1  pm,  as  observed  in  the  interstellar 
medium.  For  the  energy  of  activation,  we  will  use  the  value  by 
Ramseier  [1967]  of  59.8  kJ  mol-1,  following  the  approach 
of  Goldsby  and  Kohlstedt  [2001], 

6.4.  Input  to  the  Andrade  Model 

[95]  The  Andrade  model  (52)  contains  two  dimensionless 
empirical  parameters,  a  and  [3.  The  parameter  (3  quantifies 
the  density  and  mobility  of  the  defects  determining  the 
anelasticity  of  the  material.  As  seen  from  (53),  the  parameter 
has  units  s  “  Pa  *.  Emergence  of  a  parameter  with  fractional 
dimensions  is  an  inconvenience,  and  one  should  presume  that 
(3  is  a  fractional  power  of  a  product  or  ratio  of  parameters  with 
less  exotic  dimensions.  The  experimental  literature  contains 
few  constraints  on  the  value  of  / 3  in  general.  No  such  data  for 
ices  has  been  published  so  far.  This  leaves  us  no  other  way  but 
to  rely  on  the  data  published  for  silicates.  This  line  of  rea¬ 
soning  is  legitimate,  taken  the  remarkably  universal  nature  of 
the  Andrade  model  and  its  applicability  to  so  many  solids.  In 
Figure  4,  for  each  experiment  depicted,  we  have  represented 
(3  as  a  function  of  the  viscosity  ?y  and  shear  modulus  p.  There, 
the  fits  obtained  on  microcreep  data  are  distinguished  from 
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Andrade  Parameter  (3 


0.5x1  O'11  1.0x1  O'11  1.5X10'11 


Figure  4.  Attenuation  properties  of  olivine.  The  Andrade  parameter  p  is  depicted  as  a  function  of  the 
product  fj7{  1  ~~  a)rfa,  where  a  is  the  Andrade  exponential,  /i  is  rigidity,  and  77  is  the  viscosity.  These  data 
are  based  on  spectra  inferred  from  forced  oscillation  (diamonds)  and  microcreep  (triangles)  data  obtained 
by  Tan  et  al.  [2001]  (black)  and  by  Jackson  et  al.  [2002]  (grey).  Fits  to  the  data  by  a  least  square  method 
yield  the  black  and  grey  curves  corresponding  to  the  Tan  et  al.  and  the  Jackson  et  al.  data,  respectively. 
Results  are  discussed  in  the  text. 


those  inferred  from  spectra  of  forced  oscillations.  The  two 
methods  are  complementary,  in  that  the  frequency  bands  of 
their  maximal  sensitivity  are  different.  Besides,  analysis  of 
microcreep  data  offer  higher  resolution  when  the  viscosity 
decreases.  We  fitted  each  series  of  experiments,  using  the 
method  of  least  squares.  Overall,  the  two  data  sets  are  con¬ 
sistent  with  each  other.  Both  are  fitted  well  by  a  linear  curve 
whose  slope  is  close  to  1.  Specifically,  the  slope  is  ~1.02 
for  the  data  set  from  Jackson  et  al.  [2002]  and  is  —  1 .03  for 
the  data  set  from  Tan  et  al.  [2001].  The  R-squared  value 
is  0.62  for  the  former  set  but  is  only  0.39  for  the  latter  set. 
The  difference  is  largely  due  to  the  greater  scattering  of  the 
microcreep  data. 

[96]  Over  the  considered  range  of  frequencies,  the  vis¬ 
cosity  decreases  by  about  two  orders  of  magnitude,  while  the 
parameter  (3  decreases  by  one  order  of  magnitude.  Despite 
some  scattering  in  the  values  of  P,  we  can  observe  the  fol¬ 
lowing  relationship  between  this  parameter  and  the  param¬ 
eters  and  rj: 


times,  such  as  the  extended  Burgers  element  [e.g.,  McCarthy 
and  Castillo -Rogez,  2011]  and  the  model  by  Cole  [1995].  We 
cannot  rule  out  the  possibility  that  that  relationship  involves 
an  extra  factor,  which  would  have  been  close  to  one,  though. 

[97]  At  this  stage,  this  is  the  only  insight  on  the  quantifi¬ 
cation  of  (3  that  appears  available.  It  should  certainly  be 
employed  with  caution  outside  the  range  of  conditions  at 
which  the  aforementioned  data  sets  had  been  acquired  (that  is 
for  p  in  the  range  1 0  1  °—  1 0  1 3  s-Q  Pa-1).  We  thus  adopt 
equation  (63)  as  a  first-order  approximation,  keeping  in  mind 
that  that  relationship  may  have  to  be  updated  when  a  richer 
volume  of  laboratory  data  becomes  available.  However,  for 
the  data  at  hand,  relationship  (63)  is  valid  within  a  factor  of 
two.  We  would  emphasize  that  relationship  (63)  is  only  meant 
for  diffusion-creep-driven  anelasticity.  It  does  not  necessar¬ 
ily  apply  to  other  mechanisms,  especially  to  those  involving 
nonlinear  strain. 

[98]  With  (63)  taken  into  consideration,  the  complex 
compliance  can  be  expressed  through  the  Maxwell  time: 


=  AT(1-a)  n-a  =  J1 


(63) 


Ax) 


J[  1  +  (*7 mx)  01  r  (1  +  <*)] 


i 

VX 


(64) 


The  interrelation  between  the  viscosity  and  P  is  not  surprising 
at  all.  Both  depend  on  the  density  of  defects  and  their 
mobility,  and  on  the  appropriate  activation  energy  [e.g.,  Webb 
and  Jackson,  2003].  It  should  also  be  recalled  that  a  depen¬ 
dency  of  the  form  (tu/)“  with  a  <  1  is  typical  of  other  models 
that  aim  to  represent  the  complex  superposition  of  relaxation 


In  the  case  of  diffusion  creep,  the  Maxwell  time  rM  =  ripi  can 
be  inferred  from  equation  (60)  as: 


tm 


3RGTmd2 

42F„,p 


(65) 
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where  we  neglect  Dh,  as  the  creep  is  volume  diffusion  only. 
Pioneered  by  Raj  [1975],  expression  (65)  found  its  confir¬ 
mation  in  experiments  on  grain  boundary  diffusion  in  rocks 
[Tan  et  ai,  2001;  Jackson  et  al,  2002;  Webb  and  Jackson, 
2003].  Expressions  (63)  and  (65)  indicate  that  the  parame¬ 
ter  ff  is  sensitive,  mainly  through  viscosity,  to  the  temperature 
and  the  grain  size.  If  the  values  of  the  stresses  caused  by 
steady  state  convection  and  by  tide  differ  by  several  orders 
of  magnitude  (which  is  believed  to  be  the  case  of  Enceladus 
and  Europa),  then  different  mechanisms  may  be  involved  in 
accommodating  these  stresses  (for  example,  volume  diffu¬ 
sion  and  grain  boundary  sliding,  respectively).  In  such  a  sit¬ 
uation,  convection  and  dissipation  may  still  be  linked  to  the 
material  defect  geometry  (such  as  grain  size),  though  to  a 
lesser  extent,  since  grain  boundary  sliding  depends  on  the 
grain  size  with  an  exponent  of  only  1.4  [Goldsbv  and 
Kohlstedt,  2001]. 

[99]  For  the  Andrade  parameter  a,  a  range  of  values 
between  1/5  and  1/2  has  been  reported  most  frequently  in  the 
literature  on  attenuation  in  ice  [McCarthv  and  Castillo  -Rogez, 
2011], 

[100]  In  the  general  case  of  self-diffusion  in  solids,  theo¬ 
retical  study  by  Lifshits  and  Shikin  [1965]  suggest  that  a 
should  be  close  to  0.5  in  the  case  of  grain  boundary  diffusion. 
Jellinek  and  Brill  [1956],  too,  report  the  value  of  0.5,  though 
in  their  paper  there  is  not  enough  information  to  confirm  that 
the  transient  response  was  indeed  driven  by  the  volume  dif¬ 
fusion  creep.  Lee  and  Morris  [2010]  investigated  the  micro¬ 
physical  reasons  determining  the  value  of  a  for  grain 
boundary  diffusion,  and  established  that  a  is  primarily  a 
function  of  the  presence  of  impurities  at  the  grain  boundaries 
as  well  as  irregularities  in  the  grain  shape.  As  volume  diffu¬ 
sion  is  also  a  function  of  the  grain  size,  we  suspect  that 
nonuniform  grain  geometry  plays  a  role  in  increasing  the 
value  of  a.  Besides  the  above  references,  we  have  little 
constraint  on  the  dissipation  driven  by  the  volume  diffu¬ 
sion.  Thus,  we  explore  a  wide  range  of  values  for  a,  from 
0.2  to  0.5. 

[101]  The  above  treatment  addressed  an  anelastic  regime,  in 
which  diffusion  of  defects  is  reversible,  so  deformation  is 
recoverable.  In  the  viscosity-driven  regime,  though,  diffusion 
is  accompanied  by  boundary  sliding,  which  is  not  recover¬ 
able.  While  the  attenuation  rate  in  the  latter  regime  is 
inversely  proportional  to  the  viscosity,  Raj  and  Ashby  [1971] 
argued  that  description  of  this  process  by  the  Maxwell  model 
would  require  some  adjustment  of  the  parameters.  As  their 
theory  was  developed  for  a  certain  grain  geometry,  the 
applicability  of  their  results  to  icy  material  remains  to  be 
confirmed  by  further  theoretical  and  experimental  research. 

6.5.  Attenuation  Modeling  Approach:  Summary 

[102]  Within  the  combined  attenuation  model,  we  take  into 
account  proton  reorientation,  viscosity,  and  anelasticity  of  the 
material.  While  no  interaction  between  proton  reorientation 
and  viscosity  or  anelasticity  has  been  reported,  each  mecha¬ 
nism  dominates  the  attenuation  response  over  a  specific  range 
of  temperatures,  so  the  overlap  of  the  mechanisms  at  any 
given  time  is  minimal  anyway.  Thus,  we  assume  the  two 
mechanisms  to  act  in  parallel.  To  this  end,  we  obtain  the 
overall  compliance  function  by  summing  up  the  compliances 
appropriate  to  each  of  the  mechanisms  involved. 


[103]  In  the  case  of  proton  reorientation,  we  infer  the 
components  of  the  compliance,  input  to  the  numerical  inte¬ 
gration  code,  through 


Tm[J{x)\ 

Ke[Ax)] 


tan  S 


2  [tan/)] 


TP  X 

1  +  Tp  X2  ' 


(66) 


The  relaxation  time  i>  is  furnished  by  equation  (58),  while 
[tanb]max  is  calculated  with  aid  of  (59),  using  parameters 
presented  in  Table  1.  Most  of  these  being  based  on  the 
experimental  work  by  Tatibouet  et  al.  [1981],  we  thus  have 
only  one  data  set.  As  we  also  lack  constraints  on  the  effec¬ 
tive  modulus  of  the  material  in  that  regime,  we  assume  that 
Re(J)  ~  J.  Further  research  on  this  topic  is  necessary  in  order 
to  confirm  the  measurements  of  Tatibouet  et  al.  [1981].  In  the 
case  under  consideration,  the  relative  contribution  from  that 
mechanism  to  the  overall  despinning  evolution  is  limited,  so 
that  uncertainty  in  the  values  of  the  parameters  showing  up 
in  (58)  bears  little  implication  on  the  overall  results.  Still,  we 
included  this  mechanism  for  the  sake  of  completeness,  in 
order  to  highlight  a  process  that  could  bear  greater  implica¬ 
tions  in  other  contexts. 

[104]  The  combined  contribution  from  the  viscous,  elastic, 
and  anelastic  attenuation  is  computed  within  a  composite 
approach,  which  assumes  a  Maxwell  behavior  at  frequencies 
below  the  threshold  defined  by  (56)  and  an  anelastic-driven 
regime  at  higher  frequencies,  using  the  Andrade  model,  and 
rendered  by  (53)  and  (54).  However,  it  should  be  men¬ 
tioned  that  while  equation  (56)  suggests  a  simple  relationship 
between  Xo  and  Tm,  in  reality  the  transitional  region  from 
the  anelastic  regime  to  the  viscosity-dominated  regime  is 
complex  [Jackson  et  al,  2002]  and  may  encompass  several 
frequency  decades. 

[105]  The  temperature  enters  equation  (54)  through  the 
temperature  dependence  of  the  viscosity,  which  is  computed 
with  aid  of  equation  (60).  The  unrelaxed  shear  modulus  /i, 
too,  depends  upon  the  temperature,  but  this  dependence  is 
weak.  So  we  set  /i  =  3.3  GPa  [Parameswaran,  1987].  Other 
relevant  parameters  are  listed  in  Table  1. 

[106]  These  parameters  are  established  for  pure  water  ice, 
which  is  an  assumption  generally  adopted  for  icy  satellites, 
which,  like  lapetus,  are  dominated  by  water  ice.  Still,  the 
presence  of  soluble  (e.g.,  ammonia)  and  insoluble  (rock) 
impurities  are  expected  to  affect  the  attenuation  behavior  of 
the  material. 

[107]  Their  impact  depends  on  their  distribution  at  the  grain 
boundaries  (see  McCarthy  and  Castillo -Rogez  [2011]  for  a 
review).  Small  polar  molecules  may  also  be  incorporated  in 
the  lattice.  Dopants,  such  as  ammonia,  are  generally  known  to 
increase  dissipation,  and  to  shift  the  proton-reorientation- 
caused  dissipation  peaks  toward  a  lower  temperature  [e.g., 
Petrenko  and  Whitworth,  1999;  Oguro,  2001],  but  we  lack 
experimental  data  in  conditions  relevant  to  lapetus,  which  are 
needed  to  properly  quantify  that  effect. 

[108]  Soluble  impurities  tend  to  increase  the  grain  boundary 
width  w  and  to  create  regions  whose  thermodynamic  prop¬ 
erties  depart  from  those  of  a  pure  water  ice.  As  a  result, 
impurities  can  create  a  local  region  whose  viscosity  is  lower 
than  that  of  the  water  ice  grains.  In  this  case,  the  assumption 
that  boundary  diffusion  is  not  contributing  to  the  material 
creep,  as  generally  assumed  in  geophysical  modeling  does 
not  apply. 
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Figure  5.  (top)  The  dependence  of  the  material  viscosity  and  of  the  relaxation  time  r  upon  the  temperature 
and  the  grain  size  d  in  the  case  of  diffusion  creep  relevant  to  this  study  (see  equation  (60)).  The  shaded  region 
corresponds  to  the  realistic  range  of  values  assumed  by  the  period  of  tidal  forcing,  in  the  course  of  lapetus’ 
despinning.  The  arrow  indicates  the  approximate  viscosity  at  which  the  convection  onset  is  expected 
[Robuchon  etal.,  2010].  A  purely  viscoelastic  dissipation  model  implies  that  at  high  viscosities  the  material 
responds  elastically,  while  at  low  values  of  viscosities  the  viscosity  is  the  dominating  factor  driving  the 
response,  (bottom)  The  dependence  of  the  parameter  /?  on  the  temperature,  grain  size,  and  a,  computed  from 
equation  (63).  The  unrelaxed  shear  modulus  is  taken  as  3.3  GPa  \Parameswaran,  1987].  (The  parameter  /3 
has  units  s~a  Pa  '.) 


[109]  Although  we  currently  lack  information  systematic 
enough  for  modeling  the  impact  of  soluble  impurities  on 
material  attenuation  at  the  frequencies  relevant  to  lapetus’ 
history,  we  suspect  that  the  effect  may  play  a  major  role  in 
assisting  grain  boundary  sliding  by  boundary  diffusion  or 
boundary  migration.  We  shall  simulate  that  effect  by  testing 
different  values  for  the  activation  energy  of  defects  at  the 
grain  boundaries. 

6.6.  Numerical  Illustration  on  Simple  Models  of  lapetus 

[110]  We  present  several  numerical  runs,  to  illustrate  how 
various  processes  described  in  section  6.3  could  have  con¬ 
tributed  to  the  behavior  of  a  homogeneous  water-ice  lapetus. 


[111]  The  dependence  of  tm  on  temperature,  calculated 
through  fonnulae  (65)  and  (61),  is  plotted  in  Figure  5  (top). 
In  the  case  of  lapetus,  the  value  of  tm  becomes  of  the  same 
order  as  the  tidal  forcing  period  after  formation,  when  the 
viscosity  becomes  greater  than  ~1014  Pa  s.  As  pointed  out  by 
Castillo -Rogez  et  al.  [2007]  and  Robuchon  et  al.  [2010], 
convection  is  likely  to  begin  while  the  material  viscosity  is 
at  least  one  order  of  magnitude  greater  than  the  said  value, 
a  circumstance  which  slows  down  the  tidal  despinning. 
Flowever,  once  the  rotation  period  is  of  the  order  of  100  h,  a 
convective  lapetus  gets  more  dissipative  and  more  prone  to 
faster  tidal  evolution.  It  then  becomes  a  crucial  step  to  explain 
how  lapetus’  rotation  period  increased  from  about  10  h  to 
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Figure  6.  A  toy  model  of  dissipation  in  a  “simplified  lapetus”  consisting  of  a  pure  water  ice  and  having  a 
homogeneous  temperature  distribution,  computed  from  equation  (54).  (a)  The  inverse  of  sin  e  as  a  function 
of  the  temperature  and  a,  for  a  grain  size  of  100  //m  and  an  initial  spin  period  of  10  h.  (b)  The  same  as 
Figure  6a  except  for  an  initial  spin  period  of  100  h.  White  curves  indicate  the  viscosity  (in  Pa  s,  logarithm 
value),  (c)  The  same  as  Figure  6a  except  assuming  an  activation  energy  for  the  grain  boundary  defect 
diffusion  to  be  40  kJ  mol1,  instead  of  the  value  of  49  kJ  mol1  accepted  by  Goldsby  and  Kohlstedt  [2001]. 
The  arrows  point  to  a  viscosity  of  2  x  1 0 1 5  Pa  s  at  about  which  convection  onset  is  expected  [after  Robuchon 
et  al.,  2010].  This  example  reveals  the  importance  of  the  presence  of  second-phase  volatile  impurities 
decreasing  the  grain  boundary  viscosity. 


about  100  h,  whereafter  the  tidal  dissipation  could,  in  theory, 
be  sufficient  to  drive  despinning  to  completion. 

[112]  By  (63),  we  express  (3  as  a  function  of  a,  which  is  the 
other  parameter  entering  the  Andrade  model  (52).  In  Figure  5 
(bottom),  we  depict  the  resulting  ranges  of  values,  for  dif¬ 
ferent  temperatures,  assuming  the  material  to  be  a  pure  water 
ice.  Figure  5  (bottom)  tells  us  that  in  the  case  of  water  ice,  /3 
takes  values  lower  than  ~  10~10  s~“  Pa-1,  for  the  interval  of 


temperatures  that  interior  of  lapetus  is  expected  to  have  had 
over  its  evolution. 

[113]  Based  on  Figure  5  (bottom),  we  chart  the  value  of 
(sin|e|)— 1  for  lapetus  (using  the  approach  described  in  section 
3.5  and  equation  (54)),  assuming  the  satellite  is  homogenous 
and  composed  of  pure  water  ice,  for  a  spin  period  of  10  h 
and  100  h  (Figures  6a  and  6b).  When  a  tends  toward  zero  the 
behavior  of  the  material  tends  toward  a  Maxwell  body. 
Taking  a  equal  to  zero  is  equivalent  to  assuming  that  the 
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Figure  7.  Dissipation  in  the  “simplified  lapetus”  consisting  of  a  pure  water  ice  and  having  a  homogeneous 
temperature  distribution.  The  color  background  shows  the  dissipation  rate  characteristic  k2  sin  e  as  a  function 
of  the  grain  size  and  of  the  temperature.  The  contours  correspond  to  the  logarithm  of  the  material  viscosity, 
in  Pa  s.  The  calculation  has  been  made  for  a  spin  period  of  10  h  and  a  =  0.3.  White  curves  indicate  the  vis¬ 
cosity  (in  Pa  s,  logarithm  value).  The  grey  region  corresponds  to  conditions  where  the  model  employed  in 
this  study  becomes  unphysical  due  to  the  presence  of  melt. 


geometry  and  distribution  of  the  defects  accommodating  the 
tidal  stress  are  uniform  in  the  material,  a  situation  that  is 
obviously  not  expected  in  nature.  As  a  increases,  the 
broadening  of  the  relaxation  peaks  in  the  material  results  in 
increasing  the  temperature  range  at  which  dissipation  can  be 
significant.  The  feature  observed  in  the  bottom  right  comer 
of  Figures  6a  and  6b,  for  temperatures  lower  than  180  K,  is 
the  signature  of  proton  reorientation,  which  provides  a  lower 
bound  to  (sin|e|)  1  but  is  significant  only  at  short  periods.  On 
these  plots,  arrows  indicate  the  regime  of  temperatures  for 
which  convection  onset  is  expected,  after  Robuchon  et  al. 
[2010],  The  main  criterion  adopted  in  this  study  is  that  the 
viscosity  becomes  smaller  than  2  x  1015  Pa  s.  This  implies  a 
bound  on  the  maximal  temperature  reached  in  the  material. 
The  bound  varies  as  a  function  of  grain  size,  from  1 80  K  for 
d  =  1 0~6  m  to  230  K  for  d  =  1 0~4  m.  Within  the  Maxwell  body 
approach  (i.e.,  for  (3  — »  0  ),  the  temperature  has  to  approach 
the  water-ice  melting  temperature,  for  dissipation  to  be  large 
enough  to  promote  despinning  [cf.  Castillo -Rogez  et  al., 
2007].  On  the  other  hand,  in  the  present  modeling,  the 
value  of  1/sin  e  can  become  as  small  as  ~10  at  the  time  of 
convection  onset,  especially  for  values  of  a  greater  than  0.2, 
which  may  be  sufficient  to  drive  despinning. 

[114]  The  dependence  of  dissipation  on  the  grain  size  is 
further  illustrated  by  Figure  7.  Figure  7  shows  the  value  of 
k2  sin  e  as  a  function  of  the  grain  size  and  temperature.  The 
product  k2  sin  e  determines  the  rate  of  tidal  damping.  Figure  7 
demonstrates  that  dissipation  increases  with  the  decrease  of 
the  grain  size. 

7.  Modeling  and  Results 

7.1.  Putting  Everything  Together:  The  Overall 
Architecture  of  the  Model 

[115]  In  the  current  study,  we  focused  on  lapetus’  tidal 
response,  while  the  geophysical  part  of  the  model  largely 


repeated  the  one  developed  by  Castillo -Rogez  et  al.  [2007], 
with  some  further  amendments  regarding  the  parameters  of 
the  short-lived  radioisotope  decay,  as  suggested  by  Castillo  - 
Rogez  et  al.  [2009].  The  latter  study  suggested  a  range  of 
times  of  formation  of  3.4  to  5.4  Myr  after  the  production  of 
calcium  aluminum  inclusions  (CAIs),  so  that  heating  from 
26A1  decay  results  in  compacting  most  of  the  interior  porosity, 
a  condition  sine  qua  non  for  preservation  of  the  large  depar¬ 
ture  of  lapetus’  shape  from  hydrostatic  equilibrium.  In  the 
present  study  we  assumed  a  time  of  formation  of  ~4  Myr  after 
CAIs,  leading  to  a  scenario  consistent  with  the  current  non¬ 
hydrostatic  shape  of  lapetus.  Originally,  26A1  came  about  as  a 
part  of  the  silicate  grains  that  accreted  into  the  planetesimals 
that  later  formed  icy  satellites.  The  concentration  of  26A1 
being  subject  to  exponential  decay,  the  time  after  CAIs 
defines  how  much  26A1  is  left  by  the  beginning  of  lapetus’ 
accretion  (which  is  a  fast  process  believed  to  have  lasted  for 
only  about  105  years).  Thus,  the  time  after  CAIs  determines 
the  initial  concentration  of  26A1  considered  in  the  model. 
For  details,  see  Castillo -Rogez  et  al.  [2009]. 

[1 16]  In  our  study,  we  assumed  the  initial  spin  period  to 
lie  between  7  and  1 1  h.  As  demonstrated  by  Castillo -Rogez 
et  al.  [2009],  lapetus’  semimajor  axis  did  not  evolve,  over 
the  history  of  the  satellite,  by  more  than  a  few  thousand 
kilometers.  Hence,  in  the  current  study  we  neglected  its 
variations.  The  moment  of  inertia  C  was  recomputed  at  each 
time  step  as  a  function  of  the  appropriate  values  of  the  density 
profile  and  shape.  The  necessity  to  take  the  evolution  of  C 
into  account  stems  from  the  fact  that  the  equatorial  and  polar 
radii  could  have  evolved  by  more  than  15%  over  the  time 
span  between  lapetus’  formation  and  freezing  of  its  shape 
[Castillo -Rogez  et  al.,  2009]. 

[117]  It  would  be  important  to  emphasize  that  the  term 
“fossilization”  implies  stabilization  of  the  overall  geometrical 
shape  of  the  satellite.  Fossilization  took  place  when  the  spin 
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Figure  8.  A  “typical”  scenario  of  evolution  of  Iapetus’  geophysical  and  dynamical  properties,  computed 
under  the  following  set  of  assumptions:  the  time  of  accretion  is  4  Myr  after  the  production  of  CAls,  the  initial 
spin  period  is  1 1  h,  the  value  of  the  Andrade  parameter  a  is  0.3 .  The  evolution  of  six  parameters  is  furnished: 
(d)  the  tidal  phase  lag  e  =  e22oo  (corresponding  to  the  principal  tidal  frequency  x  —  X2200X  the  inverse  sine  of 
the  lag,  and  the  Love  number  k2,  (c)  the  maximum  temperature  achieved  in  the  model;  (b)  the  ratio  between 
the  spin  frequency  u>  and  the  mean  orbital  motion  n;  and  (a)  the  equatorial  radius. 
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Figure  9.  The  despinning  timescale  for  lapetus,  as  a  function  of  the  initial  spin  period,  of  the  value  of  the 
Andrade  parameter  a,  and  of  various  assumptions  on  material  properties.  The  number  and  letter  in  paren¬ 
theses  next  to  each  curve  refer  to  Figures  lOa-lOd  showing  the  relaxation  time  of  ice  for  the  appropriate 
case. 


period  was  15-16  h.  This  was  when  the  lithosphere  became 
thick  enough  to  support  large  nonhydro  static  anomalies.  This 
process  did  not  involve  freezing  of  the  whole  interior.  While 
the  thermal  properties  before  and  after  fossilization  were 
different,  fossilization  was  not  a  special  landmark  from  the 
viewpoint  of  the  applicability  of  our  geophysical  assump¬ 
tions.  Specifically,  the  Andrade  model  was  equally  applicable 
before  and  after  fossilization.  We  in  fact  assume  this  model  to 
work  since  the  end  of  the  accretion,  when  the  spin  period  was 
7-11  h.  We  also  assume  that  the  virtual  lack  of  porosity 
excludes  additional  dissipation  in  the  brittle  regime  (cracks, 
friction  between  fragments). 

[l is]  From  Castillo -Rogez  et  al.  [2009],  we  borrowed  a 
module  computing  the  complex  tidal  Love  number  k2  for 
multilayered  objects,  following  the  method  described  by 
Castillo  et  al.  [2000]  and  Tobie  et  al.  [2005].  This  module 
yields  the  time  evolution  of  the  following  quantities:  the 
material  phase  lag  6  (for  each  layer),  the  Love  number,  the 
equatorial  and  polar  radii,  the  moment  of  inertia,  porosity 
compaction,  and  the  possible  differentiation.  The  knowledge 
of  the  rheology  for  each  layer  was  then  used  to  calculate  the 
overall  complex  Love  number  k2(x)  for  the  body  as  a  whole, 
and  to  find  its  dependence  upon  the  tidal  frequency  x ■  In 
this,  as  explained  in  section  3.5,  we  followed  the  standard 
procedure  by  Takeuchi  and  Saito  [1972].  Then  the  values 
of  k2  sin  e22oo(xX  for  the  current  value  of  the  tidal  frequency 
X,  were  inserted  in  fonnulae  from  section  4.2  to  simulate 
despinning. 

[119]  We  simulate  the  consequences  of  convection  by 
assuming  that  the  viscosity  never  increases  beyond  1015  Pa  s 
[Robuchon  et  al.,  2010].  In  reality,  this  constraint  is  some¬ 
what  simplistic,  since  Robuchon  et  al.  [2010]  obtained 


warmer  conditions  in  some  of  their  models.  Thus,  our  result 
yields  an  upper  bound  on  the  despinning  timescale.  In  any 
case,  the  low  dissipation  rate  resulting  from  spin  evolution  is 
not  expected  to  affect  convection  (G.  Robuchon,  personal 
communication,  2010).  We  also  took  into  account  that  the 
time  of  convection  onset  at  temperatures  greater  than  200  K  is 
of  the  order  of  several  million  years  [ Robuchon  et  al.,  2010]. 

[120]  The  parameters  used  in  the  model  are  gathered  in 
Table  1. 

7.2.  A  Typical  Scenario 

[121]  A  typical  scenario  of  despinning,  in  the  context  of 
thermal  evolution,  is  presented  in  Figure  8.  Within  this 
example,  despinning  is  achieved  over  ~1.6  Gyr,  while  the 
internal  temperature  remains  lower  than  210  K.  Over  this 
entire  time,  the  phase  lag  e  obeys  l/sin|e|  <  100.  Being  small 
after  accretion,  the  parameter  l/sin|  e|  <  100  increases  due  to 
the  low  temperature.  The  contribution  of  proton  reorientation 
at  these  low  temperatures  and  high  forcing  frequencies  keeps 
1/sinj  e|  below  100.  After  a  few  hundred  million  years  after 
formation,  we  get  l/sin|e|  ~30.  As  the  spin  period  increases, 
the  regime  of  dissipation  evolves  from  purely  anelastic  to 
viscoelastic,  and  l/sin|e|  decreases  to  about  10  and  lower,  as 
the  spin  period  becomes  greater  than  1 000  h  shortly  before  the 
end  of  despinning  (for  this  particular  scenario,  at  ~1.4  Gyr 
after  formation).  An  important  implication  is  that  despinning 
evolves  gradually,  and  the  phase  lag  evolution  is  primarily 
a  consequence  of  the  period  evolution.  This  contrasts  with 
a  Maxwell  body-based  model,  in  which  little  despinning 
occurs  until  warm  temperatures  are  achieved,  and  despinning 
occurs  on  a  very  short  timescale  (i.e.,  the  phase  lag  evolution 
is  primarily  driven  by  the  temperature  changes).  Thus,  the 
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Figure  10.  Relaxation  times  for  lapetus’  material  computed  from  equation  (65),  first  by  considering  the 
volume  diffusion  (VD)  and  grain  boundary  diffusion  (BD)  separately,  and  then  by  summing  them  (curve 
labeled  VD+BD),  after  equation  (60),  and  using  the  input  parameters  presented  in  Table  1.  The  cases  in 
Figures  1  Oa-1 0c  assume  pure  water  ice,  with  different  grain  sizes:  (a)  d  =  1 00  /zm;  (b)  d  =  1 0  //m;  and  (c)  d  = 
1  //m.  (d)  Case  assumes  the  presence  of  a  lower  activation  energy  for  impurities  located  at  the  grain 
boundaries  ( Ea  =  40  kJ  mol1).  The  shaded  regions  correspond  to  a  range  of  spin  periods  relevant  to  this 
study. 


despinning  runaway  model  and  rapid  shape  evolution  pro¬ 
posed  by  Castillo  -Rogez  et  al.  [2007]  are  not  valid.  The  latter 
model  suggested  a  change  in  the  equatorial  radius  by  almost 
30%  over  a  few  million  years.  In  the  current  study,  the  shape 
evolves  as  a  rate  almost  constant  until  its  freezing,  which 
takes  place  at  about  600  Myr  after  formation. 

7.3.  The  Despinning  Timescale 

[122]  The  dependencies  of  the  despinning  timescale  on  the 
key  input  parameters  (the  initial  spin  period,  a,  and  the  grain 
size)  are  charted  in  Figure  9.  We  consider  two  assumptions 
about  the  nature  of  the  ice  grain  boundaries.  Most  cases 
depicted  on  Figure  9  assume  that  the  activation  energy  of 
grain  boundary  defect  migration  has  the  value  suggested  by 
Goldsby  and  Kohlstedt  [2001],  e.g.,  49  kJ  mol-1.  In  some 
situations,  though,  a  lower  value  of  40  kJ  mol  1  may  be 
possible,  as  argued  by  Goodman  et  al.  [1981],  Such  a  low 
value  may  be  attributed  to  the  presence  of  salt  impurities.  It  is 


beyond  the  scope  of  this  study  to  discuss  the  validity  of  the 
Goldsby  and  Kohlstedt  [2001]  versus  Goodman  et  al.  [1981] 
numbers.  The  primary  goal  of  these  numerical  runs  is  to 
explore  how  a  simulated  departure  from  the  theoretical  ther¬ 
modynamic  properties  for  pure  water  ice  affects  the  response 
of  the  material  to  mechanical  forcing. 

[123]  In  the  case  ofthe  “theoretical”  pure  water  ice,  Figure  9 
shows  that  the  despinning  time  ranges  from  about  900  to 
3700  Myr.  The  latter  is  an  upper  bound  reached  for  a  wide 
range  of  conditions.  Shorter  despinning  times  are  possible  for 
a  narrow  set  of  conditions,  namely  for  a  grain  size  lower  than 
1 0-4  m,  an  initial  spin  period  longer  than  1 0  h,  and  a  smaller 
than  0.4.  It  is  important  to  note  that  if  we  neglect  the  grain 
boundary  diffusion  creep,  and  build  the  attenuation  model 
only  using  volume  diffusion  creep,  then  the  despinning  time 
would  always  be  of  the  order  of  3.6  Gyr.  This  can  be 
understood  by  considering  Figure  10  that  shows  the  relax¬ 
ation  time  for  grain  boundary  diffusion  and  volume  diffu- 
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sion  as  a  function  of  the  temperature  and  grain  size.  For  d  > 
0.01  mm,  the  relaxation  time  of  the  material  is  determined  by 
the  volume  diffusion  for  almost  the  entire  relevant  tempera¬ 
ture  range.  As  the  grain  size  decreases,  the  grain  boundary 
diffusion  mechanism  may  dominate  the  global  response  of 
the  material  subject  to  short-period  forcing,  leading  to  a  sit¬ 
uation  where  the  volume  diffusion  drives  steady  state  creep, 
while  the  grain  boundary  diffusion  dominates  at  the  forcing 
periods  comparable  to  Iapetus’  spin  period.  If  the  grain  size  is 
of  the  order  of  0.1  mm,  the  despinning  time  shows  little 
dependence  on  the  initial  spin  period  and  a,  and  varies  by  less 
than  50  Myr  over  the  probed  parameter  range.  Decreasing  the 
grain  size  to  0.001  mm  results  in  decreasing  the  despinning 
timescale  by  almost  a  factor  4  if  the  spin  period  is  11  h. 
However,  for  a  spin  period  lower  than  9  h,  the  despinning 
timescale  is  systematically  of  the  order  of  3.6  Gyr.  In  any 
case,  with  this  approach,  we  do  not  find  despinning  times  of 
a  few  hundred  million  years  as  previously  suggested  by 
Castillo  -Rogez  etal.  [2007]  and  Robuchon  et  al.  [2010]  based 
on  rather  arbitrary  choices  on  the  input  parameters  to  the 
attenuation  model. 

[124]  The  despinning  timescale  may  be  significantly 
reduced  if  the  effect  of  impurities  is  simulated  by  decreasing 
the  energy  of  activation  and  increasing  the  grain  boundary 
width.  If  the  activation  energy  is  decreased  by  15%,  as 
illustrated  in  Figure  9,  the  despinning  timescale  may  be  less 
than  1  Gyr.  In  that  case,  the  dominant  mechanism  driving 
defect  motion  changes  from  the  volume  diffusion  to  the  grain 
boundary  diffusion,  as  the  temperature  increases  (Figure  10). 
Figure  6c  illustrates  that  for  the  same  conditions  as  repre¬ 
sented  in  Figure  6a  attenuation  is  increasing  over  a  wider 
range  of  temperatures.  For  example,  for  a  grain  size  of  1 0  //m, 
the  maximal  attenuation  is  achieved  at  a  temperature  lower  by 
20  K  than  that  for  pure  water  ice.  As  a  result,  attenuation  may 
become  geophysically  significant  for  a  range  of  temperature 
where  convection  dominates  heat  transfer,  so  that  convection 
is  not  a  constraint  on  despinning  duration  anymore.  Although 
the  choice  of  parameters  for  that  particular  test  is  arbitrary,  it 
reflects  a  widespread  configuration  that  associates  different 
mechanisms  acting  in  response  to  forcing  exerted  over  a 
broad  range  of  frequencies.  This  way,  shorter  despinning 
times  can  be  made  possible  by  exploring  a  wider  range  of 
parameters  characterizing  Iapetus’  material  chemistry. 

7.4.  Tidal  Heating 

[125]  Within  the  considered  model,  most  heat  is  produced 
deep  below  the  surface.  At  the  temperature  about  150  K 
and  a  very  low  porosity,  the  thermal  conductivity  of  an  ice- 
rock  mixture  can  be  as  large  as  3  W  (m  K.)  1  at  150  K.  As 
despinning  gets  accomplished  over  a  timescale  of  hundreds 
of  millions  of  years  (at  least),  we  conclude  that  within  the 
low-porosity  model,  tidal  friction  does  not  contribute  con¬ 
siderably  to  heating  of  the  interior,  and  thus  plays  no  role 
in  the  endogenic  activity.  (As  explained  in  Appendix  A  to 
Castillo -Rogez  et  al.  [2007],  in  the  absence  of  heat  transfer, 
the  tidal-stress-generated  temperature  increase,  over  hundreds 
of  millions  of  years,  would  be  about  15  K.) 

7.5.  The  Role  of  Iapetus’  Triaxiality 

[126]  A  spinning  satellite  is  always  subject  to  two  torques 
exerted  upon  it  by  its  host  planet.  Besides  the  aforementioned 


tidal  torque,  there  exists  a  torque  caused  by  the  satellite’s 
triaxiality.  This  torque  taken  into  account,  expression  (44)  for 
despinning  rate  should  be  amended  to: 

3  Ml ;c7k2sin|e220o|  Rl  3  B  -  A  2a3 

9  — - - - L  — - n  -T- sin 2(9  —  v) 

2  £  Mprim  a6  2  C  r3  K  ’ 

+  0(e2sine)  +  0(/2sine),  (67) 

7  standing  for  Newton’s  gravitational  constant;  0,  9,  n,  and  v 
denoting  Iapetus’  sidereal  angle,  spin  rate,  mean  motion, 
and  true  anomaly;  while  A  <  B  <  C  being  Iapetus’  principal 
moments  of  inertia  (recomputed  at  each  time  step).  Through 
its  history,  Iapetus’  stays  mostly  homogeneous  in  density, 
with  no  chemical  differentiation  and  almost  no  porosity;  so 
the  evolution  of  the  moments  of  inertia  is  defined  primarily 
by  the  change  of  shape.  After  Iapetus  gets  fossilized,  the 
moments  of  inertia  stay  constant. 

[127]  Since  for  solid  moons  and  planets  their  triaxiality 
typically  exceeds  the  height  of  the  tidal  bulge,  the  triaxiality  - 
caused  torque  exceeds  the  tidal  torque.  Moreover,  whenever 
the  orbital  eccentricity  is  not  zero,  the  triaxiality-caused 
torque  has  to  be  taken  into  account  in  modeling  the  approach 
to  the  1 : 1  spin-orbit  resonant  state,  as  was  pointed  long  ago 
by  Goldreich  [1966].  However,  during  the  despinning  stage, 
the  triaxiality-caused  torque  gets  averaged  out  and  there¬ 
fore  plays  no  major  role  in  the  rotational  dynamics,  until  the 
rotator  approaches  the  final  state.  This  is  why,  in  the  first 
approximation,  this  torque  is  often  neglected.  A  more  care¬ 
ful  study,  though,  indicates  that  sometimes  the  rotational 
behavior  becomes  more  complicated  [ Wisdom  et  al.,  1984; 
Melnikov  and  Shevchenko,  2010],  and  the  triaxiality-caused 
torque  plays  an  important  role  in  it,  especially  in  the  vicinities 
of  the  low-order  spin-orbit  resonances  like  5:2,  3:2,  etc. 

[128]  Here  we  limit  ourselves  to  several  numerical  runs  of 
equation  (67),  carried  out  by  means  of  numerical  integrator 
RA 1 5  authored  by  Everhart  [  1 9 8 5] .  In  our  opinion,  these  runs 
illustrate  the  role  of  the  triaxiality-caused  term,  at  least  in  the 
context  of  the  chosen  dissipation  model. 

[129]  We  begin  with  a  more  careful  investigation  of  the 
model  considered  by  Aleshkina  [2009].  In  that  model,  k2  sin 
£2200  was  replaced  with  k2!Q,  where  k2  was  set  constant, 
while  Q  scaled  as  the  inverse  to  the  tidal  frequency.  As 
we  already  mentioned,  this  dissipation  law  is  not  supported 
by  experimental  measurements.  The  mechanical  model  with 
k2  =  const  and  Q  oc  x  1  can  be  shortly  written  down  as  k2/Q  oc 
X ■  Branding  this  model  as  unphysical,  we  should  keep  in 
mind  that  the  more  realistic  (Andrade)  model  (53),  too,  leads 
to  a  similar  relation  k2  sin  (>2200  oc  X2200  between  the  Love 
number,  the  principal  tidal  frequency  X2200.  and  the  appro¬ 
priate  phase  lag  (>2200.  as  demonstrated  in  Appendix  A.  The 
fundamental  difference,  however,  lies  in  the  fact  that  the 
unphysical  model  extends  the  dependence  k2IQ  oc  x  to 
the  entire  range  of  frequencies  x>  including  those  corre¬ 
sponding  to  commensurabilities  other  than  1:1.  The  realistic 
model,  on  its  part,  simplifies  to  k2  sin  e22oo  oc  X2200  only  in  the 
closemost  vicinity  of  the  1 : 1  resonance  (where  X2200  — 1 ►  0), 
while  elsewhere  (including  the  vicinities  of  the  other  com¬ 
mensurabilities)  it  yields,  with  a  good  degree  of  approxima¬ 
tion:  k2  sin  e22oo  oc  X2200  (for  high  frequencies)  and  k2  sin 
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Figure  IE  Reexamination  of  the  model  employed  by  Aleshkina  [2009].  The  model  includes  both  the 
triaxiality-caused  torque  and  the  tidal  torque  (equation  (67))  and  assumes  that  the  quality  factor  entering  the 
tidal  torque  scales  as  Q  cx  x_1.  Here  the  said  model  is  implemented  by  three  numerical  runs,  with  the  same 
initial  condition  90  =  1  .In  for  the  spin  rate  but  with  slightly  different  initial  conditions  for  the  sidereal  angle: 
90=  1.5680,  1.5700,  1.5719  rad. 


£2200  oc  X2200  a)  (for  low  frequencies),  both  a  and  (1  -  a) 
being  positive. 

[130]  The  second  aspect  of  Aleshkina’s  model,  which 
we  intend  to  reexplore,  is  the  role  of  triaxiality  at  various 
initial  conditions.  Individual  evolutionary  paths,  which  differ 


through  tiny  changes  in  initial  conditions,  are  all  equally  valid 
statistical  representatives  of  the  actual  system.  So  we  reex¬ 
amined  Aleshkina’s  results  by  performing  three  simulations 
with  slightly  distinct  initial  values  of  the  sidereal  angle:  9  0  = 
1.5680, 1.5700,  1.5719  radian.  To  simplify  the  comparison  of 
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Figure  12.  Comparison  of  two  scenarios  of  Iapetus’  tidal  despinning.  One  scenario  is  based  on  the 
Andrade  dissipation  model  explored  in  this  paper;  another  is  based  on  the  model  Q  cx  x_1  employed  by 
Aleshkina  [2009].  In  both  cases,  the  triaxiality-generated  torque  was  included.  This  torque  did  not  cause 
major  changes  to  the  scenario  based  on  the  Andrade  model.  However,  it  brought  a  temporary  entrapment  (in 
the  3:4  resonance)  into  the  scenario  based  on  the  model  Q  cx  %_1.  This  happened  because  the  latter  model 
implies  an  unrealistic  increase  of  the  Q  factor  at  low  frequencies,  leading  to  an  unrealistic  decrease  of  the 
despinning  rate  and,  as  a  result,  to  a  higher  probability  of  temporary  entrapments  in  intermediate  resonances. 
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the  thus  obtained  results  with  those  reported  by  Aleshkina 
[2009],  we  began  all  three  simulations  with  the  same  initial 
spin  rate  90  =  l.7n,  and  used  for  all  the  involved  parameters 
the  same  values  as  Aleshkina  [2009].  To  decrease  the  simu¬ 
lation  time,  we  multiplied  the  tidal  torque  by  a  factor  of 
one  thousand.  This  method  is  often  employed  to  study  the 
behavior  of  dynamical  systems.  It  is  legitimate  as  long  as  the 
timescales  of  dissipative  mechanisms  much  exceed  those  of 
nondissipative  forces.  As  the  tidal  torque  is  dissipative  and 
acts  on  very  long  timescales,  we  artificially  accelerate  its 
action.  On  the  other  hand,  the  triaxiality-caused  torque  is  left 
untouched,  since  we  do  not  want  to  change  the  physics  on  the 
short  timescale. 

[131]  Figure  11  shows  the  three  resulting  scenarios  of 
lapetus’  despinning,  within  Aleshkina’s  model  implemented 
under  three  slightly  different  initial  conditions  imposed  on  9. 
The  first  choice  of  the  initial  condition,  90  =  1 .5  6  80  rad,  leads 
to  lapetus’  crossing  the  low-order  spin-orbit  resonances,  with 
no  obvious  time  delay.  The  other  two  choices,  90  =  1 .5700  rad 
and  90  =  1.57  1  9  rad,  result  in  lapetus  getting  stuck  for  a 
significant  amount  of  time  in  both  the  5:4  and  3:2  resonances. 
While  the  friction  law  Q  oc  x_1  employed  by  Aleshkina 
[2009]  guarantees  exceedingly  long  tidal-despinning  times, 
there  may  be  a  possibility  that  the  afore  described  delays  in 
low-order  resonances  led  to  further  elongation  of  those  times 
of  Aleshkina  [2009], 

[132]  To  emphasize  the  role  that  the  frequency  dependence 
of  dissipation  plays  in  the  tidal-despinning  process,  in 
Figure  12  we  compare  two  scenarios  based  on  different 
dissipation  models.  One  scenario  is  based  on  the  observation- 
based  dissipation  model  explored  in  the  present  work,  another 
on  the  model  Q  oc  employed  by  Aleshkina  [2009].  Since 
in  the  observation-based  model  the  quality  factor  does 
not  increase  rapidly  at  low  frequencies,  this  model  keeps 
despinning  fast  even  when  the  ratio  9/n  becomes  of  order 
unity.  However,  in  the  Q  oc  x  1  model  the  quality  factor 
grows  at  low  frequencies,  thus  providing  much  lower 
despinning  rate,  and  increasing  the  chances  for  getting  tem¬ 
porarily  stuck  at  the  intermediate  resonances  due  to  the  tri- 
axiality-generated  torque;  see  the  long-term  entrapment  in 
the  4:3  resonance  in  Figure  12. 

8.  Conclusions 

[133]  We  have  examined  the  despinning  history  of 
lapetus,  using  the  mechanical  model  of  Andrade  which  is 
quite  different  from  models  employed  in  the  literature  hitherto. 
Our  approach  is  based  on  the  observation  that  under  the  stress 
and  temperature  conditions  relevant  to  lapetus  following 
accretion,  the  most  likely  anelastic  and  viscoelastic  dissipation 
mechanisms  involve  the  motion  and  rearrangement  of  defects 
in  the  ice  lattice  (self-lattice  diffusion).  In  this  sense,  lapetus 
is  typical,  in  that  most  of  the  other  large  icy  moons,  too,  are 
subject  to  low  tidal  stressing.  (Internal  friction  in  satellites 
subject  to  high  tidal  stress,  such  as  Europa,  Enceladus, 
Miranda,  Ariel,  is  likely  to  be  dominated  by  dislocations  or 
grain  boundary  sliding  as  primary  drivers  of  dissipation  [e.g., 
McCarthy  and  Castillo -Rogez,  2011,  Figure  1].) 

[134]  Our  study  is  based  on  the  Andrade  model,  which  has 
proven  adequate  to  describe  anelasticity  in  a  wide  variety  of 
materials,  including  ices.  An  important  improvement  of  this 


approach  upon  earlier  models  is  that  it  has  enabled  us  to 
combine  the  viscous  and  transient  creep  components  within 
one  model,  in  a  self-consistent  manner.  We  also  explored  the 
role  of  proton  reorientation  as  an  additional  source  of  atten¬ 
uation  leading  to  a  possible  departure  of  the  dissipation  law 
from  the  Andrade  model  at  certain  temperatures.  Using  the 
available  empirical  data  as  well  as  necessary  extrapolations, 
we  have  demonstrated  that  tidal  dissipation  in  lapetus  is 
sufficient  to  achieve  efficient  despinning  at  temperatures 
much  below  the  water  ice  melting  temperature.  Under  these 
conditions,  the  effect  of  convection  as  a  factor  limiting  dis¬ 
sipation  and  despinning,  pointed  out  by  Castillo -Rogez  et  al. 
[2007]  and  explored  by  Robuchon  et  al.  [2010],  is  no  longer 
an  obstacle. 

[135]  At  the  same  time,  our  study  has  highlighted  certain 
difficulties  inherent  in  this  type  of  problems.  Specifically,  if 
we  model  dissipation  under  the  assumption  that  lapetus  con¬ 
sists  of  a  pure  water  ice,  the  resulting  despinning  timescales 
come  up  to  about  4  Gyr,  a  result  that  may  be  inconsistent  with 
the  available  geological  data.  Then  impurities,  whose  pres¬ 
ence  in  the  ice-dominated  satellites  is  suspected  but  whose 
affect  is  generally  not  quantified,  may  have  a  significant 
impact  on  the  attenuation  mechanism.  Rock  impurities  are 
expected  to  inhibit  boundary  sliding  at  these  stress  levels 
[Raj,  1975],  while  low-eutectic,  second-phase  volatiles  are 
likely  to  decrease  the  viscosity  locally  by  promoting  defect 
mobility.  While  the  prospect  of  investigating  these  processes 
in  future  is  most  appealing,  virtually  no  experimental  data 
are  available  so  far.  This  has  limited  our  ability  to  narrow 
down  our  estimates  for  lapetus’  despinning  timescale.  While 
the  obtained  estimate  places  the  despinning  time  within  the 
interval  from  0.9  Gyr  through  3.7  Gyr,  a  more  exact  estimate 
will  remain  unavailable  until  we  learn  more  about  the  influ¬ 
ence  of  impurities  upon  dissipation  in  ices.  (While  the 
experimental  literature  offers  some  data  on  dissipation  in 
dirty  ices,  those  data  are  typically  obtained  under  forcing 
several  orders  of  magnitude  higher  than  the  tidal  stress 
expected  in  most  icy  satellites.  Hence  the  applicability  of  such 
data  in  our  study  remains  questionable.) 

[136]  We  have  also  pointed  out  that  by  adding  the  triaxi¬ 
ality-caused  torque  to  the  tidal  one,  we  encounter  a  chaotic 
behavior  at  the  final  stage  of  despinning,  a  behavior  that 
sometimes  includes  long-term  entrapments  in  the  interme¬ 
diate  resonances.  Although  this  phenomenon  needs  further 
investigation,  our  numerical  runs  indicate  that  entrapment 
becomes  more  likely  when  one  employs  the  unphysical 
model  Q  oc  x1  >n  calculation  of  the  tidal  torque.  Such  a 
model  entails  an  unrealistic  increase  of  the  Q  factor  at  low 
frequencies,  leading  to  an  unrealistic  decrease  of  the 
despinning  rate  and,  as  a  result,  to  a  higher  probability  of 
temporary  entrapments  in  intermediate  resonances.  On  the 
other  hand,  when  calculation  of  the  tidal  torque  is  based  on  a 
more  realistic  mechanical  model,  the  tidal  torque  remains 
large  enough  even  at  the  latest  stages  of  despinning.  So  the 
despinning  rate  stays  sufficiently  high,  thus  reducing  the 
chances  of  getting  stuck  at  an  intermediate  resonance. 

[137]  Finally,  in  Appendix  A,  we  explained  why  the  dissi¬ 
pation  models  Q  ~  x“,  with  a  positive  a,  do  not  lead  to 
infinities  in  the  expressions  for  the  tidal  torque  or  force. 
(Incorrect  claims  of  existence  of  such  infinities  appeared 
in  the  literature  more  than  once.)  Although  Efroimsky  and 
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Williams  [2009]  already  addressed  this  point,  here  we 
explained  it  more  carefully. 

Appendix  A:  Behavior  of  |&2(;\;)|sin  e2(x)  = 

- Xm  [£(x)l  in  Various  Frequency  Bands, 
for  the  Andrade  and  Maxwell  Models 

[138]  Efroimskv  and  Williams  [2009]  emphasized  that  the 
empirical  model  Q  ~  x“  with  a  positive  a  does  not  doom  the 
Impq  term  of  the  torque  to  explode  on  arrival  at  the  hnpq 
resonance.  While  the  frequency  Ximpq  approaches  zero,  it 
should  be  remembered  that  the  said  tidal  term  actually  con¬ 
tains  not  ki/Q  but  A7  sin  the  phase  lag  being  related  to  the 
quality  factor  via  (6)  or  (7),  dependent  on  how  exactly  the 
quality  factor  is  defined. 

[139]  Efroimskv  and  Willianis  [2009]  mentioned  that  in  the 
limit  of  vanishing  frequency,  the  Love  number  scales  in  the 
leading  order  as 

ki(x)  ~  k2( 0)  cos  e2(x)-  (Al) 


Andrade  body,  the  complex  compliance  J  is  given  by  (53).  Its 
imaginary  and  real  parts  are 

lrn[J{x)\  =  -  —  -x_Q^sm^)r(a+  1),  (A4) 

VX  v  2  ' 

Ke[J{x)]=J  +  x~a/3cos(®P)r(a  +  1)  (A5) 

Introducing  the  Maxwell  time 

tm  =  -  =  r)J,  (A6) 

d 

and  recalling  that  A/  1,  we  can  easily  derive  from  (63) 
and  (A4)-(A5)  that 

|£;(x)|sine;(x)  ™  2^_  ^  ^sjn(f)r(a  +  >)M'“ 

for  x>Lv1>  (A7) 


However,  our  justification  of  (Al)  was  less  than  solid.  We 
referred  to  a  very  empirical  hydrodynamical  treatment  by 
Alexander  [1973].  We  also  mentioned,  with  no  proof,  that  by 
using  formulae  from  Churkin  [1998]  it  is  possible  to  prove 
(Al)  for  a  broad  class  of  models. 

[140]  Here  we  shall  demonstrate  that  (Al)  works  well  for 
the  Maxwell  model.  For  the  Andrade  body,  though,  it  gets 
replaced  with  a  more  complicated  expression.  Nonetheless, 
as  proven  below,  crossing  of  a  resonance  portends  no  diffi¬ 
culties,  because  the  crucial  factor  k 7  sin  e/  vanishes  smoothly 
as  the  moon  gets  into  a  resonance. 

[141]  Generally,  |L/x)|sin  efx)  =  -Xm[F/(x)].  It  is  straight¬ 
forward  from  (27)  that 


™>-WT)Jb5TAJ 

3 _ E.e[J(x)\  +  ilm[J(x)\ 

2(l-l)lZe[J(x)}  +  A,J  +  il  m[J(X)] 


so  that 


ki(x)  = 


2(1-1) 


{Ke[J(x)\f  +  (Im[J(x)\y  +  A,JTle[J{X)}  +  iA,  J  Im[J(X)\ 
(1Ze[J(x)\  +  A,  J)2  +  (lm[J(X)]f 

(A2b) 

wherefrom 


\ki(x)\ sine,(x)  =  -lm[k,(x)\  = 


2(1-1) 

_ -A,  Jlm[J{x)\ _ 

(TZe[J(x)\  +  A,jf  +  (lm[J{X)))2' 


(A3) 


where  J=  J(0)  =  1/p  —  Vp(0)  may  be  chosen  as  the  unrelaxed 
compliance  (see  the  paragraph  after  formula  (26)).  For  an 


3  1 

l*/(x)|sme,(x)  ~  2(/ -  1)7/^™^' 

for  tm'  >  x  >  rMlAf,  (A8) 

3 

|fc/(x)|sine/(x)  ~  ^A,tmX  for  rM14[‘»x.  (A9) 

Expressions  (A7)  and  (A8)  coincide  with  the  frequency 
dependencies  of  |J(x)|sin  <5(x)  =  ~lm[J(x)\  in  the  high-  and 
low-frequency  bands,  as  can  be  demonstrated  from  (53). 
However,  (A9)  renders  a  feature  characteristic  of  the  tidal 
lagging,  and  not  of  that  in  the  material.  We  owe  this  feature  to 
self-gravitation,  i.e.,  to  the  presence  of  the  first  term  in  the 
denominator  of  formula  (25). 

[142]  For  realistic  values  of  the  parameters  of  icy  satellites 
(say,  p  ~  109  Pa,  r/~  101 1  Pa  s,  and  ~  103),  the  condition 
t~m  Ay1  x  puls  x  below  10~5  Hz.  This  indicates  that  \kfx)\ 
sin  e/(x)  behaves  as  (A9)  only  in  an  extremely  close  vicinity 
of  the  resonance  corresponding  to  vanishing  of  the  mode  X- 
Still,  it  is  important  that  the  quantity  |£/(\)jsin  <y(x)  grows  to 
a  finite  maximum  and  then  vanishes  rapidly  but  smoothly,  as 
the  frequency  falls  to  zero.  So  the  expression  for  the  tidal 
torque  or  tidal  force  experiences  no  infinities  when  the  moon 
approaches  a  resonance. 

[143]  As  mentioned  in  section  3.3,  in  reality  the  potential  U 
and  therefore  also  k  are  functions  not  of  the  positively  defined 
frequencies  \  but  of  the  tidal  modes  u>  which  can  be  positive 
or  negative.  Employing  \  instead  of  ui,  we  must  compensate 
for  this  abuse  by  multiplying  the  lag  efximpq),  “by  hand”, 
with  sgno Jimpq.  This  way,  in  actual  computations  each 
expression  (A7)-(A9)  must  be  amended  with  this  factor. 
Then,  for  example,  (A9)  will  read  as  \kix,mpq)\sin  efximpq )  ~ 
777  [7  ^  C MXimpq  sgn tO/mpq  477  TT  A  C From  this  we 
see  that  most  naturally,  the  Impq  term  of  the  torque  changes 
its  sign  on  crossing  the  Impq  resonance. 

[144]  Finally,  we  would  point  out  that  the  general  expres¬ 
sion  (A2a)  furnishes  the  expression  for  |^/(x)|>  while  (A4) 
and  (A5)  render  the  expression  for  tan  e/x)  =  -Tm[ki(X)\/lZe 
TZelkfxf],  From  these  expressions,  it  is  easy  to  demonstrate 
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that  (Al)  indeed  works,  for  the  Maxwell  model,  in  the  low- 
frequency  limit.  However^  for  the  Andrade  model  the  rela¬ 
tion  between  |k/(x)|  and  |k/(0)|  becomes  more  complicated. 
Interestingly,  for  the  Kelvin- Voigt  model  relation  (Al)  is 
exact  at  all  frequencies. 

Notation 

a  semimajor  axis. 

A,  B,  C  moments  of  inertia  of  the  satellite  (A  <  B  <  C). 

Co  initial  concentration  in  Bjerrum  defects  at  a  refer¬ 
ence  temperature. 
d  ice  grain  size. 

Db  diffusivity  of  grain  boundary  diffusion. 

Dv  diffusivity  of  volume  diffusion. 
e  eccentricity  of  the  secondary. 

E  tidal  energy. 

Ea  activation  energy  for  diffusion  creep. 

Ef  activation  energy  of  Bjerrum  defect  formation. 

Em  activation  energy  of  Bjerrum  defect  migration. 
Flmp(i)  inclination  functions. 

g  surface  gravity. 

Gipq{e)  eccentricity  polynomials. 

i  inclination  of  the  secondary. 

I  degree  (spherical  harmonics). 

J,  .7(0)  compliance. 

kt  tidal  Love  number  of  degree  I. 

Mprim  mass  of  the  primary  (satellite). 

Msec  mass  of  the  secondary  (planet). 

M  mean  anomaly. 
n  mean  motion. 

p  slope  of  the  attenuation  spectrum. 

Pi  Legendre  polynomial  of  degree  /. 

Q  dissipation  factor. 
r  radius. 

R  satellite’s  mean  radius. 

Rg  gas  constant. 
t  time. 

T  temperature. 

Tm  melting  temperature. 
u  shear  strain. 

Vm  molar  volume  of  defects. 
w  grain  boundary  thickness. 

W  stationary  tidal  change  of  the  potential. 
a,  (3  input  parameters  to  the  Andrade  creep  equation. 

X  physical  frequencies  of  deformation. 

6  material  phase  lag. 

At  time  lag. 

A E  energy  loss. 

A  dissipation  strength, 
e  tidal  phase  lag. 
r]  steady  state  viscosity. 

7  Newton’s  gravitational  constant. 

A  longitude. 

r  gamma  function. 

/i,  ^(0)  unrelaxed  shear  modulus. 
v  true  anomaly. 
lu  argument  of  the  pericenter. 

longitude  of  the  node. 
p  initial  phase. 

</>  latitude. 
a  shear  stress. 


tm  Maxwell  time. 
tv  Voigt  time. 

tp  relaxation  time  for  the  proton  reorientation 
mechanism. 

T  tidal  torque. 

0  Heaviside  function. 

9  sidereal  angle. 

9  spin  rate. 

£  dimensionless  mean  moment  of  inertia. 

(,  v  tensor  indices. 
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